Method and apparatus for calculating blood flow parameters

ABSTRACT

A method for determining tissue type includes quantitatively determining a tissue blood flow (TBF) by deconvoluting Q(t) and C a (t), where Q(t) represents a curve of specific mass of contrast, and C a (t) represents an arterial curve of contrast concentration, and quantitatively determining a tissue blood volume (TBV) by deconvoluting Q(t) and C a (t). The method also includes quantitatively determining a tissue mean transit time (TMTT) by deconvoluting Q(t) and C a (t), and quantitatively determining a tissue capillary permeability surface area product (TPS) by deconvoluting Q(t) and C a (t). The method also includes determining a tissue type based on the TBF, the TBV, the TMTT, and the TPS.

CROSS REFERENCE TO RELATED APPLICATIONS

[0001] This application claims the benefit of U.S. provisionalapplication No. 60/243,196 filed Oct. 25, 2000.

BACKGROUND OF THE INVENTION

[0002] This invention relates generally to methods and apparatus forimaging systems and, more particularly, for calculating blood flowparameters.

[0003] Stroke is a leading cause of disability among adults in NorthAmerica. A National Institute for Neurological Diseases and Stroke(NINDS) trial was the first study to demonstrate that thrombolytictreatment within the first 3 hours of symptom onset is beneficial tostroke victims. However, treating patients beyond the initial 3 hours ofillness with a thrombolytic treatment may carry a risk of intracranialbleeding which has resulted in a setback in the use of thrombolysis.This setback in the use of thrombolysis has prompted the realizationthat a more rigorous selection of patients according to certain criteriamay reduce the risk of intracranial bleeding and thus, extend thetherapeutic window for this therapy beyond the current 3-hour limit tobenefit more patients.

BRIEF SUMMARY OF THE INVENTION

[0004] In one aspect, a method for determining tissue type is provided.The method includes quantitatively determining a tissue blood flow (TBF)by deconvoluting Q(t) and C_(a)(t), where Q(t) represents a curve ofspecific mass of contrast, and C_(a)(t) represents an arterial curve ofcontrast concentration, and quantitatively determining a tissue bloodvolume (TBV) by deconvoluting Q(t) and C_(a)(t). The method alsoincludes quantitatively determining a tissue mean transit time (TMTT) bydeconvoluting Q(t) and C_(a)(t), and quantitatively determining a tissuecapillary permeability surface area product (TPS) by deconvoluting Q(t)and C_(a)(t). The method also includes determining a tissue type basedon the TBF, the TBV, the TMTT, and the TPS.

[0005] In another aspect, a system including at least one of a computedtomography system and a nuclear magnetic resonance system is provided.The system is configured to quantitatively determine a tissue blood flow(TBF) by deconvoluting Q(t) and C_(a)(t), where Q(t) represents a curveof specific mass of contrast, and C_(a)(t) represents an arterial curveof contrast concentration, and quantitatively determine a tissue bloodvolume (TBV) by deconvoluting Q(t) and C_(a)(t). The system is alsoconfigured to quantitatively determine a tissue mean transit time (TMTT)by deconvoluting Q(t) and C_(a)(t), and quantitatively determine atissue capillary permeability surface area product (TPS) bydeconvoluting Q(t) and C_(a)(t). The system is also configured todetermine a tissue type based on the TBF, the TBV, the TMTT, and theTPS.

[0006] In another embodiment, a computer readable medium encoded with aprogram executable by a computer for processing scanning data isprovided. The program is configured to instruct the computer toquantitatively determine a tissue blood flow (TBF) by deconvoluting Q(t)and C_(a)(t), where Q(t) represents a curve of specific mass ofcontrast, and C_(a)(t) represents an arterial curve of contrastconcentration, and quantitatively determine a tissue blood volume (TBV)by deconvoluting Q(t) and C_(a)(t). The program is also configured toinstruct the computer to quantitatively determine a tissue mean transittime (TMTT) by deconvoluting Q(t) and C_(a)(t), and quantitativelydetermine a tissue capillary permeability surface area product (TPS) bydeconvoluting Q(t) and C_(a)(t). The program is also configured toinstruct the computer to determine a tissue type based on the TBF, theTBV, the TMTT, and the TPS.

BRIEF DESCRIPTION OF THE DRAWINGS

[0007]FIG. 1 is a pictorial view of a CT imaging system.

[0008]FIG. 2 is a block schematic diagram of the system illustrated inFIG. 1.

DETAILED DESCRIPTION OF THE INVENTION

[0009] At present, there are both experimental and clinical data tosuggest that the level of residual blood flow in the ischemic zone couldbe a useful indicator of the risk of intracranial bleeding.

[0010] The importance of absolute measurement in stroke is alsopredicated by the fact that there exist distinct thresholds of cerebralblood flow (CBF) for various functions of the brain. It is possible touse these thresholds to decide whether a particular brain region issalvageable, hence thrombolytics should be used to restore CBF, or isalready infarcted such that thrombolysis would achieve little butinstead would increase the risk of intracranial bleeding.

[0011] The brain has an intricate system of control to maintain cerebralnormal limits when cerebral perfusion pressure decreases but remainswithin certain limits. This is achieved by dilation of the resistanceblood vessels, i.e. the arterioles, vascular resistance and an increaseof the cerebral blood volume (CBV). The concept of autoregulation of CBFcan be used to advantage in stroke. For ischemic but viable tissue,autoregulation would lead to increased CBV so that mean transit time(MTT) is prolonged since MTT is the ratio of CBV and CBF because of theCentral Volume Principle. On the other hand for ischemic but nonviabletissue, autoregulation is abolished such that both CBV and CBF arereduced but MTT may remain normal or only slightly elevated.

[0012] In summary, the absolute measurements of CBF, CBV and MTT wouldpermit the distinction between salvageable and infarcted tissue by thefollowing criteria: Ischemic Tissue Type CBF CBV MTT Viable −− + ++Infarcted −− − −/+

[0013] It is the mismatch between CBF and CBV that discriminatessalvageable and infarcted tissue. The corollary of this is thatmeasurement of CBF alone will not reliably differentiate between viableand non-viable ischemic tissue. In addition, the ability to monitorquantitative changes in CBV may obviate the need to performsupplementary tests to evaluate the ‘vascular’ reserve of a strokepatient. In brief, the classical test of ‘vascular’ reserve evaluatesthe change in CBF in response to increases in pCO2 (partial pressure ofCO2) in brain tissue which can be induced either by increasing inspiredair concentration of CO2 or by intravenous administration of a drug,such as a carbonic anhydrase inhibitor such as, for example, Diamox. Fora normal brain, raised pCO2 in tissue would induce a large increase inCBF. For ischemic tissue, where autoregulation is still intact, becauseCBV is already activated, the increase in CBF will be attenuated. Thus,tissue with a positive reserve is viable whereas that with little or noreserve left is at risk of infarction.

[0014] As used herein Cerebral Blood Flow (CBF) is the volume flow ofblood through vasculature including the large conductance vessels,arteries, arterioles, capillaries, venules, veins and sinuses. Ittypically has units of ml/min/100 g. Also as used herein, Cerebral BloodVolume (CBV) is the volume of blood in the vasculature including thelarge conductance vessels, arteries, arterioles, capillaries, venules,veins and sinuses. It typically has units of ml/g. Additionally, as usedherein Mean Transit Time (MTT) references that blood traverses thevasculature through different pathlengths such that there does not existan unique transit time from the arterial inlet to the venous outlet.Instead there is a distribution of transit times and the mean transittime is the mean time of such a distribution. Furthermore, MinimumTransit Time (TTmin), as used herein, is the minimum time intervalbetween entry at the arterial inlet and exit at the venous outlet ofblood or contrast medium. The Central Volume Principle relates the abovethree quantities in the following relationship:${CBF} = {\frac{CBV}{MTT}.}$

[0015] Referring to FIGS. 1 and 2, a multi-slice scanning imagingsystem, for example, computed tomography (CT) imaging system 10, isshown as including a gantry 12 representative of a “third generation” CTimaging system. Gantry 12 has an x-ray source 14 that projects a beam ofx-rays 16 toward a detector array 18 on the opposite side of gantry 12.Detector array 18 is formed by a plurality of detector rows (not shown)including a plurality of detector elements 20 which together sense theprojected x-rays that pass through an object, such as a medical patient22. Each detector element 20 produces an electrical signal thatrepresents the intensity of an impinging x-ray beam and hence theattenuation of the beam as it passes through object or patient 22.During a scan to acquire x-ray projection data, gantry 12 and thecomponents mounted thereon rotate about a center of rotation 24. FIG. 2shows only a single row of detector elements 20 (i.e., a detector row).However, a multislice detector array 18 includes a plurality of paralleldetector rows of detector elements 20 so that projection datacorresponding to a plurality of parallel slices are, or can be acquiredsimultaneously during a scan.

[0016] Rotation of gantry 12 and the operation of x-ray source 14 aregoverned by a control mechanism 26 of CT system 10. Control mechanism 26includes an x-ray controller 28 that provides power and timing signalsto x-ray source 14 and a gantry motor controller 30 that controls therotational speed and position of gantry 12. A data acquisition system(DAS) 32 in control mechanism 26 samples analog data from detectorelements 20 and converts the data to digital signals for subsequentprocessing. An image reconstructor 34 receives sampled and digitizedx-ray data from DAS 32 and performs high-speed image reconstruction. Thereconstructed image is applied as an input to a computer 36 which storesthe image in a mass storage device 38.

[0017] Computer 36 also receives commands and scanning parameters froman operator via console 40 that has a keyboard. An associated cathoderay tube display 42 allows the operator to observe the reconstructedimage and other data from computer 36. The operator supplied commandsand parameters are used by computer 36 to provide control signals andinformation to DAS 32, x-ray controller 28 and gantry motor controller30. In addition, computer 36 operates a table motor controller 44 whichcontrols a motorized table 46 to position patient 22 in gantry 12.Particularly, table 46 moves portions of patient 22 through gantryopening 48. In one embodiment, computer 36 includes a device 50, forexample, a floppy disk drive or CD-ROM drive, for reading instructionsand/or data from a computer-readable medium 52, such as a floppy disk orCD-ROM. In another embodiment, computer 36 executes instructions storedin firmware (not shown). Computer 36 is programmed to perform functionsdescribed herein, accordingly, as used herein, the term computer is notlimited to just those integrated circuits referred to in the art ascomputers, but broadly refers to computers, processors,microcontrollers, microcomputers, programmable logic controllers,application specific integrated circuits, and other programmablecircuits.

[0018] In one embodiment, system 10 is used to perform CT scans anddetermines blood flow parameters such as Tissue Blood Flow (TBF), TissueBlood Volume (TBV), Tissue Mean Transit Time (TMTT) and Tissue CapillaryPermeability Surface Area Product (TPS) as described below. In anotherembodiment, a Nuclear Magnetic Resonance (NMR) system (not shown) scansand determines TBF, TBV, TMTT and TPS blood flow parameters such asdescribed below. In an exemplary embodiment, system 10 scans anddetermines cerebral blood flow parameters such as Cerebral Blood Flow(CBF), Cerebral Blood Volume (CBV), and Cerebral Mean Transit Time(CMTT).

[0019] In one embodiment, system 10 is used to determine tissue type.More specifically, a tissue blood flow (TBF) is quantitativelydetermined by deconvoluting Q(t) and C_(a)(t), where Q(t) represents thetissue residue function and is a curve of specific mass of contrast intissue, and C_(a)(t) represents an arterial curve of contrastconcentration. Also a tissue blood volume (TBV), a tissue mean transittime (TMTT), and a tissue capillary permeability surface area product(TPS) are quantitatively determined by deconvoluting Q(t) and C_(a)(t).In another embodiment, a Nuclear Magnetic Resonance (NMR) system (notshown) scans and determines tissue type by quantitatively determiningTBF, TBV, TMTT and TPS.

[0020] In one embodiment, the arterial curve of contrast concentrationmeasured by system 10 and NMR system (not shown) is corrected forpartial volume averaging as described herein. For example, during acranial scan, but not limited to cranial scans, arterial regions withinthe vascular territories of the cerebral arteries (anterior and middle)are identified, and used to generate the measured arterial curve ofcontrast concentration, C_(a)(t). The measured arterial curve ofcontrast concentration is related to the arterial curve of contrastconcentration, C_(a)(t) by C_(a)′(t)=k*C_(a)(t), where k is the partialvolume averaging scaling factor as explained in greater detail below. Avenous region either within the sagittal or tranverse sinuses islocated, and C_(v)(t) is generated where C_(v)(t) C_(a)(t)*h(t) and h(t)is the transit time spectrum of the brain, as explained herein.C_(a)′(t) and C_(v)(t) are deconvolved to find $\frac{h(t)}{k}.$

[0021] And a trailing slope of C_(a)′(t) is extrapolated with amonoexponential function to find C_(a,ex)′(t) which is convolved with$\frac{h(t)}{k}$

[0022] to find C_(v,ex)(t). Where k is a partial volume averaging (PVA)scaling factor and is determined according to$k = {\frac{\int_{0}^{\infty}{{C_{a,{ex}}^{\prime}(t)}\quad {t}}}{\int_{0}^{\infty}{{C_{v,{ex}}(t)}\quad {t}}}.}$

[0023] The measured arterial curve of contrast concentration, C_(a)′(t),is then corrected for partial volume averaging by dividing with thefactor k to arrive at the arterial curve of contrast concentration,C_(a)(t).

[0024] A CBF and a CBV functional maps are generated with C_(a)′(t) anddivided by the PVA scaling factor, k.

[0025] In one embodiment, a convolution integral for a tissue residuefunction Q(t) is: Q(t) = ∫₀^(t)C_(a)(t − u)h(u)  u

[0026] where C_(a)(t) is an arterial curve of contrast concentrationcurve and h(t) is an impulse residue function. Additionally, Q(t) isequal to the convolution of C_(a)(t) and h(t) when the system is linearand stationary. Note that when C_(a)(t)=δ(t−u), then the tissue residuefunction Q(t) is: Q(t) = ∫₀^(t)δ(t − u)h(u)  u = h(t)

[0027] This relationship lends itself to the following interpretation ofthe convolution integral, e.g.${{C_{a}(t)} = {\sum\limits_{n}{{C_{a}\left( {t - {n\quad \Delta \quad t}} \right)}{\delta \left( {t - {n\quad \Delta \quad t}} \right)}}}},$

[0028] therefore the tissue residue function Q(t) is:$\quad \begin{matrix}{{{Q(t)} = {\int_{0}^{t}{{C_{a}\left( {t - u} \right)}{h(u)}\quad {u}}}},} & {{which}\quad {is}} \\{{= {\int_{0}^{t}{\sum\limits_{n}{\left\lbrack {{C_{a}\left( {t - u} \right)}{\delta \left( {t - u - {n\quad \Delta \quad t}} \right)}} \right\rbrack {h(u)}\quad {u}}}}},} & {{which}\quad {is}} \\{{= {\sum\limits_{n}{\int_{0}^{t}{{C_{a}\left( {t - u} \right)}{\delta \left( {t - u - {n\quad \Delta \quad t}} \right)}{h(u)}\quad {u}}}}},} & {{which}\quad {is}} \\{= {\sum\limits_{n}{{C_{a}\left( {n\quad \Delta \quad t} \right)}{{h\left( {t - {n\quad \Delta \quad t}} \right)}.}}}} & \quad\end{matrix}\quad$

[0029] A discretization of the tissue residue function convolutionintegral is determined by dividing a time interval [0,t] into m equalintervals of length Δt, each using the rectangular rule for quadratureof the convolution integral, according to:${Q(t)} = {{Q\left( {m\quad \Delta \quad t} \right)} = {\sum\limits_{i = 0}^{m - 1}\quad {{C_{a}\left( {\left\lbrack {m - i} \right\rbrack \Delta \quad t} \right)}{h\left( {i\quad \Delta \quad t} \right)}\Delta \quad {t.}}}}$

[0030] In a matrix notation, the convolution integral for the tissueresidue function after discretization is Q=Ah and therefore Q=Ah is:$\begin{pmatrix}{Q\left( {\Delta \quad t} \right)} \\{Q\left( {2\Delta \quad t} \right)} \\\vdots \\{Q\left( {\left( {M - 1} \right)\Delta \quad t} \right)} \\{Q\left( {M\quad \Delta \quad t} \right)}\end{pmatrix} = {\quad{\begin{bmatrix}{C_{a}\left( {\Delta \quad t} \right)} & 0 & 0 & \cdots & \cdots & 0 \\{C_{a}\left( {2\Delta \quad t} \right)} & {C_{a}\left( {\Delta \quad t} \right)} & 0 & 0 & \cdots & 0 \\\vdots & \vdots & \vdots & \vdots & \vdots & \vdots \\{C_{a}\left( {\left( {M - 1} \right)\Delta \quad t} \right)} & {C_{a}\left( {\left( {M - 2} \right)\Delta \quad t} \right)} & \cdots & \cdots & \cdots & \left. {{C_{a}\left( {M - N} \right)}\Delta \quad t} \right) \\{C_{a}\left( {M\quad \Delta \quad t} \right)} & {C_{a}\left( {\left( {M - 1} \right)\Delta \quad t} \right)} & {C_{a}\left( {\left( {M - 2} \right)\Delta \quad t} \right)} & \cdots & \cdots & {{C_{a}\left( \left( {M - N + 1} \right) \right)}\Delta \quad t}\end{bmatrix}\begin{pmatrix}{h(0)} \\{h\left( {\Delta \quad t} \right)} \\{h\left( {2\Delta \quad t} \right)} \\\vdots \\{h\left( {\left( {N - 1} \right)\Delta \quad t} \right)}\end{pmatrix}}}$

[0031] where Q is a M×1 vector, A is a M×N matrix, and h is a N×1vector. A least squares solution, ĥ, of Q=Ah facilitates minimizing thenorm of a tissue residue vector, r=Q=Aĥ according to: ∥r∥²=∥Q−Aĥ∥².

[0032] In one embodiment, an equality constraint is incorporated intothe least squares problem such that when h is the impulse residuefunction of a linear and stationary flow system, it will satisfy thetime causality, i.e. if the injection of tracer does not occur at timezero, then some beginning elements of h should be set to zero (timecausality), and defining a minimum transit time, i.e. h will have aplateau of a duration equal to the minimum transit time. The requirementfor time causality and minimum transit time can be written as thefollowing equalities: $\begin{matrix}{h_{1} = 0} \\{h_{2} = 0} \\\vdots \\{h_{I} = 0} \\{h_{I + 1} = h_{I + 2}} \\{h_{I + 2} = h_{I + 3}} \\\vdots \\{h_{I + L - 1} = h_{I + L}}\end{matrix}$

[0033] where the first I elements of h are zero, and the duration of theplateau is (L−1) dt, where dt is the sampling interval of each elementof h. The equality constraints can be written compactly as the matrixequation Ch=b: $C = {\begin{matrix}I \\\quad \\\quad \\\quad \\{L - 1} \\\quad\end{matrix}\begin{bmatrix}\overset{\overset{I}{}}{\begin{matrix}1 & 0 & \cdots & 0 \\0 & 1 & \cdots & 0 \\\vdots & \vdots & \cdots & \vdots \\0 & 0 & \cdots & 1 \\0 & 0 & \vdots & 0 \\0 & 0 & \vdots & 0 \\\vdots & \vdots & \vdots & \vdots \\0 & 0 & \cdots & 0 \\\quad & \quad & \quad & \quad\end{matrix}} & \overset{\overset{L}{}}{\begin{matrix}0 & 0 & 0 & \cdots & 0 & 0 \\0 & 0 & 0 & \cdots & 0 & 0 \\\vdots & \vdots & \vdots & \cdots & \vdots & \vdots \\0 & 0 & 0 & \cdots & 0 & 0 \\1 & {- 1} & 0 & \cdots & 0 & 0 \\0 & 1 & {- 1} & \cdots & 0 & 0 \\\vdots & \vdots & \vdots & \cdots & \vdots & \vdots \\0 & 0 & 0 & \cdots & 1 & {- 1} \\\quad & \quad & \quad & \quad & \quad & \quad\end{matrix}} & \overset{\overset{N - I - L}{}}{\begin{matrix}0 & \cdots & \cdots & \cdots & \cdots & \cdots & \cdots & \cdots & \cdots & 0 \\0 & \cdots & \cdots & \cdots & \cdots & \cdots & \cdots & \cdots & \cdots & 0 \\\vdots & \cdots & \cdots & \cdots & \cdots & \cdots & \cdots & \cdots & \cdots & \vdots \\0 & \cdots & \cdots & \cdots & \cdots & \cdots & \cdots & \cdots & \cdots & 0 \\0 & \cdots & \cdots & \cdots & \cdots & \cdots & \cdots & \cdots & \cdots & 0 \\0 & \cdots & \cdots & \cdots & \cdots & \cdots & \cdots & \cdots & \cdots & 0 \\\vdots & \cdots & \cdots & \cdots & \cdots & \cdots & \cdots & \cdots & \cdots & \vdots \\0 & \cdots & \cdots & \cdots & \cdots & \cdots & \cdots & \cdots & \cdots & 0 \\\quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad\end{matrix}}\end{bmatrix}}$

[0034] where C is a M₁×N matrix, h is the impulse residue fuifction, andb=0 is a N×1 zero vector. Note that M₁−I+L−1.

[0035] Also C can be partitioned as: $\begin{matrix}I \\\quad \\\quad \\\quad \\{L - 1} \\\quad\end{matrix}\begin{bmatrix}\overset{\overset{C_{1}}{}}{\underset{\underset{{I + L} = {1 = M_{1}}}{}}{\begin{matrix}1 & 0 & \cdots & 0 & 0 & 0 & 0 & \cdots & 0 \\0 & 1 & \cdots & 0 & 0 & 0 & 0 & \cdots & 0 \\\vdots & \vdots & \cdots & \vdots & \vdots & \vdots & \vdots & \cdots & 0 \\0 & 0 & \cdots & 1 & 0 & 0 & 0 & \cdots & 0 \\0 & 0 & \vdots & 0 & 1 & {- 1} & 0 & \cdots & 0 \\0 & 0 & \vdots & 0 & 0 & 1 & {- 1} & \cdots & 0 \\\vdots & \vdots & \vdots & \vdots & \vdots & \vdots & \vdots & \cdots & \vdots \\0 & 0 & \cdots & 0 & 0 & 0 & 0 & \cdots & 1\end{matrix}}} & \overset{\overset{C_{2}}{}}{\underset{\underset{{N - {({I + L - 1})}} = {N - M_{1}}}{}}{\begin{matrix}0 & 0 & \cdots & \cdots & \cdots & \cdots & \cdots & \cdots & \cdots & \cdots & 0 \\0 & 0 & \cdots & \cdots & \cdots & \cdots & \cdots & \cdots & \cdots & \cdots & 0 \\\vdots & \vdots & \cdots & \cdots & \cdots & \cdots & \cdots & \cdots & \cdots & \cdots & \vdots \\0 & 0 & \cdots & \cdots & \cdots & \cdots & \cdots & \cdots & \cdots & \cdots & 0 \\0 & 0 & \cdots & \cdots & \cdots & \cdots & \cdots & \cdots & \cdots & \cdots & 0 \\0 & 0 & \cdots & \cdots & \cdots & \cdots & \cdots & \cdots & \cdots & \cdots & 0 \\\vdots & \vdots & \cdots & \cdots & \cdots & \cdots & \cdots & \cdots & \cdots & \cdots & \vdots \\{- 1} & 0 & \cdots & \cdots & \cdots & \cdots & \cdots & \cdots & \cdots & \cdots & 0\end{matrix}}}\end{bmatrix}$

[0036] where C₁ is a M₁×M₁ square matrix with full rank, and C₂ is aM₁×(N−M₁) matrix with only one nonzero element, wherein the nonzeroelement (M₁+1, l)th element is equal to −1. Similarly h can bepartitioned as: ${h = \begin{pmatrix}h_{1} \\h_{2}\end{pmatrix}},$

[0037] where h₁ is a M₁×1 vector of the first M₁ elements of h₁ and h₂is a N−M₁×1 vector of the remaining elements.

[0038] With these partitions of C and h, the constraint equation, Ch=b,can be written as:${\left\lbrack {C_{1}\quad C_{2}} \right\rbrack \begin{pmatrix}h_{1} \\h_{2}\end{pmatrix}} = {\left. b\Rightarrow{C_{1}h_{1}} \right. = {\left( {b - {C_{2}h_{2}}} \right).}}$

[0039] Since C₁ is square and of full rank, its inverse exists, then theconstraint equation can be written as: h₁=C₁ ⁻¹(b−C₂h₂), where theinverse of C₁ is: $\begin{matrix}I \\\quad \\\quad \\{C_{1}^{- 1} =} \\\quad \\{L - 1}\end{matrix}\begin{bmatrix}\begin{matrix}\underset{}{I} \\\begin{matrix}1 & 0 & \cdots & 0 \\0 & 1 & \cdots & 0 \\\vdots & \vdots & ⋰ & \vdots \\0 & 0 & \cdots & 1\end{matrix}\end{matrix} & \begin{matrix}\underset{}{L - 1} \\\begin{matrix}0 & 0 & 0 & 0 & 0 \\0 & 0 & 0 & 0 & 0 \\\vdots & \vdots & \vdots & \vdots & \vdots \\0 & 0 & 0 & 0 & 0\end{matrix}\end{matrix} \\\begin{matrix}0 & 0 & 0 & 0 \\0 & 0 & 0 & 0 \\\vdots & \vdots & \vdots & \vdots \\0 & 0 & 0 & 0 \\0 & 0 & 0 & 0\end{matrix} & \begin{matrix}1 & 1 & \cdots & 1 & 1 \\0 & 1 & \cdots & 1 & 1 \\\vdots & \vdots & ⋰ & \vdots & \vdots \\0 & 0 & 0 & 1 & 1 \\0 & 0 & 0 & 0 & 1\end{matrix}\end{bmatrix}$

[0040] As a check, h₁=C₁ ⁻¹(b−C₂h₂) is evaluated to ensure that it doesnot violate the constraint described herein. For example, since b is azero vector, h₁ can be simplified as: h₁=−C₁ ⁻¹C₂h₂. As statedpreviously herein, C₂ is a M₁×N−M₁ matrix of the form:$C_{2} = {\begin{bmatrix}0 & 0 & \cdots & 0 \\0 & 0 & \cdots & 0 \\\vdots & \vdots & \vdots & \vdots \\0 & 0 & \cdots & 0 \\0 & 0 & \cdots & 0 \\0 & 0 & \cdots & 0 \\\vdots & \vdots & \vdots & \vdots \\0 & 0 & \cdots & 0 \\{- 1} & 0 & \cdots & 0\end{bmatrix} = {\left\lbrack {\begin{matrix}\overset{1}{} \\0 \\0 \\\vdots \\0 \\{- 1}\end{matrix}\begin{matrix}\underset{}{N - M_{1} - 1} \\\begin{matrix}0 & 0 & \cdots & 0 \\0 & 0 & \cdots & 0 \\\vdots & \vdots & \vdots & \vdots \\0 & 0 & \cdots & 0 \\0 & 0 & \cdots & 0\end{matrix}\end{matrix}} \right\rbrack = \begin{bmatrix}0 & 0 \\{- 1} & 0\end{bmatrix}}}$

[0041] and h₂ can be partitioned as: $h_{2} = {\begin{pmatrix}h_{M_{1} + 1} \\h_{M_{1} + 2} \\\vdots \\h_{N}\end{pmatrix} = \begin{pmatrix}h_{M_{1} + 1} \\h_{2}\end{pmatrix}}$

[0042] then, $\begin{matrix}{{C_{2}h_{2}} = {{\,_{1}^{M_{i} - 1}\begin{bmatrix}\overset{\overset{1}{}}{0} & \overset{\overset{N - M_{i} - 1}{}}{0} \\{- 1} & 0\end{bmatrix}}\begin{pmatrix}h_{M_{i} + 1} \\h_{2}\end{pmatrix}\begin{matrix}1 \\{N - M_{i} - 1}\end{matrix}}} \\{{= {\begin{pmatrix}0 \\{- h_{M_{i} + 1}}\end{pmatrix}\begin{matrix}{M_{i} - 1} \\1\end{matrix}\quad {and}}},} \\{h_{1} = {\begin{pmatrix}h_{1} \\h_{2} \\\vdots \\h_{I} \\h_{I + 1} \\h_{I + 2} \\\vdots \\h_{I + L}\end{pmatrix} = {- {C_{1}^{- 1}\begin{pmatrix}O \\{- h_{M_{i} + 1}}\end{pmatrix}}}}} \\{= {\,_{L - 1}^{1}\left\lbrack {\frac{\overset{\overset{M_{i} - 1}{}}{X}}{X}\left. \begin{matrix}\begin{matrix}\overset{\overset{1}{}}{0} \\0 \\\vdots \\0 \\\quad \\{- 1} \\{- 1} \\\vdots\end{matrix} \\{- 1}\end{matrix} \right\rbrack \begin{pmatrix}O \\{- h_{M_{i} + 1}}\end{pmatrix}\begin{matrix}{M_{i} - 1} \\1\end{matrix}} \right.}} \\{= {\left( \frac{\begin{matrix}0 \\0 \\\vdots \\0\end{matrix}}{\begin{matrix}h_{M_{i} + 1} \\h_{M_{i} + 1} \\\vdots \\h_{M_{i} + 1}\end{matrix}} \right)\begin{matrix}I \\{L - 1}\end{matrix}}}\end{matrix}$

[0043] as the time causality and minimum transit time constraintsrequired.

[0044] In one embodiment, the equality constraint Ch=b affect on thesolution of the least squares problem O=Ah can be determined accordingto: First, A is partitioned in the same way as C:${A =_{M}\left\lbrack \begin{matrix}{\overset{\overset{M_{1}}{}}{A_{1}},} & \overset{\overset{N - M_{1}}{}}{A_{2}}\end{matrix}\quad \right\rbrack},$

[0045] then Ah can be written as:${Ah} = {{\left\lbrack {A_{1},A_{2}} \right\rbrack \begin{pmatrix}h_{1} \\h_{2}\end{pmatrix}} = {{{A_{1}h_{1}} + {A_{2}h_{2}}} = {{{A_{1}\left( {{- C_{1}^{- 1}}C_{2}h_{2}} \right)} + {A_{2}h_{2}}} = {\left( {A_{2} - {A_{1}C_{1}^{- 1}C_{2}}} \right){h_{2}.}}}}}$

[0046] Therefore with the equality constraint, Ch=b, the least squaresproblem Ah=Q is equivalent to the ‘reduced’ problem of (A₂−A₁C₁⁻¹C₂)h₂=Q or in short A_(r)h₂=Q where A_(r)=(A₂A₁C₁ ⁻¹C₂).

[0047] From the least squares solution, ĥ₂, of the reduced problemA_(r)h₂=Q, the full impulse residue function, h, can be reconstitutedby: $h = \begin{pmatrix}h_{1} \\{\hat{h}}_{2}\end{pmatrix}$

[0048] and h₁ is a (I+L−1)×1 vector whose first I elements are zero andthe subsequent L−1 elements are equal to the first element of ĥ₂.

[0049] For example, for the evaluation of A₂−A₁C₁ ⁻¹C₂ ${\begin{matrix}1 \\\quad \\\quad \\{{C_{1}^{- 1}C_{2}} = {M_{1} - 1}} \\\quad\end{matrix}\begin{bmatrix}\begin{matrix}\overset{M_{1} - 1}{} \\X\end{matrix} & \begin{matrix}\overset{1}{} \\0\end{matrix} \\X & \begin{matrix}0 \\\vdots \\0 \\1 \\\vdots \\1\end{matrix}\end{bmatrix}}{\begin{matrix}\quad \\\quad \\\quad \\{M_{1} - 1} \\\quad \\\quad \\\quad \\1\end{matrix}\begin{bmatrix}\begin{matrix}\overset{1}{} \\0 \\0 \\\vdots \\\vdots \\\vdots \\0\end{matrix} & \begin{matrix}\overset{N - M_{1} - 1}{} \\\quad \\\quad \\0 \\\quad \\\quad \\\quad\end{matrix} \\{- 1} & {0\quad \cdots \quad 0}\end{bmatrix}}$ ${\begin{matrix}1 \\\quad \\\quad \\{= {M_{1} - 1}} \\\quad\end{matrix}\begin{bmatrix}\begin{matrix}\overset{1}{} \\0\end{matrix} & \begin{matrix}\overset{N - M_{1} - 1}{} \\{0\quad \cdots \quad 0}\end{matrix} \\\begin{matrix}0 \\\vdots \\0 \\{- 1} \\\vdots \\{- 1}\end{matrix} & \begin{matrix}\quad \\\quad \\0 \\\quad \\\quad \\\quad\end{matrix}\end{bmatrix}} = {\begin{matrix}\quad \\{M_{1} - L + 1} \\\quad \\\quad \\\quad \\\quad \\{L - 1} \\\quad\end{matrix}\begin{bmatrix}\begin{matrix}1 \\0 \\0 \\\vdots \\0\end{matrix} & \begin{matrix}\overset{N - M_{1} - 1}{} \\\quad \\{\begin{matrix}0 \\\quad\end{matrix}\quad} \\\quad\end{matrix} \\{- 1} & \quad \\\vdots & 0 \\{- 1} & \quad\end{bmatrix}}$

[0050] where X is a plurality of partitions of matrices that are notused. A₁ is then partitioned according to: $A_{1} = {\begin{matrix}{M - L + 1} \\{L - 1}\end{matrix}\left\lbrack \quad \begin{matrix}\overset{\overset{M_{1 - L + 1}}{}}{X} & \overset{\overset{L - 1}{}}{T_{1}} \\X & T_{2}\end{matrix}\quad \right\rbrack}$

[0051] such that:${A_{1}C_{1}^{- 1}C_{2}} = {{{\,_{L - 1}^{M - L + 1}\begin{bmatrix}\begin{matrix}\overset{M_{1} - L + 1}{} \\X\end{matrix} & \begin{matrix}\overset{L - 1}{} \\T_{1}\end{matrix} \\X & T_{2}\end{bmatrix}}{\,_{L - 1}^{M_{1} - L + 1}\begin{bmatrix}\begin{matrix}\begin{matrix}\begin{matrix}\begin{matrix}\overset{1}{} \\0\end{matrix} \\0\end{matrix} \\\vdots\end{matrix} \\0\end{matrix} & \begin{matrix}\begin{matrix}\begin{matrix}\overset{N - M_{1} - 1}{} \\\quad\end{matrix} \\0\end{matrix} \\\quad\end{matrix} \\\begin{matrix}\begin{matrix}{- 1} \\\vdots\end{matrix} \\{- 1}\end{matrix} & 0\end{bmatrix}}} = {\,_{L - 1}^{M - L + 1}\begin{bmatrix}\begin{matrix}\begin{matrix}\overset{1}{} \\\quad\end{matrix} \\{T_{1}\begin{pmatrix}{- 1} \\\vdots \\{- 1}\end{pmatrix}}\end{matrix} & \overset{\overset{N - M_{1} - 1}{}}{0} \\{T_{2}\begin{pmatrix}{- 1} \\\vdots \\{- 1}\end{pmatrix}} & 0\end{bmatrix}}}$

[0052] therefore, A₂−A₁C₁ ⁻¹C₂ is the same as A₂ except that its firstcolumn is modified by adding the M×1 vector$\,_{L - 1}^{M - L + 1}\begin{bmatrix}{T_{1}\begin{pmatrix}1 \\\vdots \\1\end{pmatrix}} \\{T_{2}\begin{pmatrix}1 \\\vdots \\1\end{pmatrix}}\end{bmatrix}$

[0053] to it.

[0054] In one embodiment, a smoothness constraint by Lagrange Multiplieris incorporated into the least squares solution of the reduced problemA_(r)h₂=Q. In fact, if h is the impulse residue function of aphysiologically realizable flow system, then its elements would changesmoothly. In one embodiment, a matrix γF is appended to A_(r)=(A₂−A₁C₁⁻¹C₂) and vector γd to Q to the reduced problem to obtain:$\begin{pmatrix}q \\{\gamma \quad d}\end{pmatrix} = {\left. {\begin{bmatrix}A_{r} \\{\gamma \quad S}\end{bmatrix}h_{2}}\Leftrightarrow p \right. = {{{Eh}_{2}\quad {where}\quad p} = \begin{pmatrix}q \\{\gamma \quad d}\end{pmatrix}}}$

[0055] is a M₂×1 vector and $E = \begin{bmatrix}A_{r} \\{\gamma \quad S}\end{bmatrix}$

[0056] is a M₂×(N−M₁) matrix, and p=Eh₂ is referred to the appendedreduced problem.

[0057] The least squares solution, ĥ₂, of the appended reduced problemp=Eh₂ facilitates minimizing the normal of the residual vector accordingto: $\begin{matrix}{{r}^{2} = \quad {\left( {\begin{pmatrix}Q \\{\gamma \quad d}\end{pmatrix} - {\begin{bmatrix}A_{r} \\{\gamma \quad S}\end{bmatrix}{\hat{h}}_{2}}} \right)^{T}\left( {\begin{pmatrix}Q \\{\gamma \quad d}\end{pmatrix} - {\begin{bmatrix}A_{r} \\{\gamma \quad S}\end{bmatrix}{\hat{h}}_{2}}} \right)}} \\{= \quad {\left( {\left\lbrack {Q^{T}\quad \gamma \quad d^{T}} \right\rbrack - {{\hat{h}}_{2}^{T}\left\lbrack {A_{r}^{T}\quad \gamma \quad S^{T}} \right\rbrack}} \right)\left( {\begin{pmatrix}Q \\{\gamma \quad d}\end{pmatrix} - {\begin{bmatrix}A_{r} \\{\gamma \quad S}\end{bmatrix}{\hat{h}}_{2}}} \right)}} \\{= \quad {{\left\lbrack {Q^{T}\quad \gamma \quad d^{T}} \right\rbrack \begin{pmatrix}Q \\{\gamma \quad d}\end{pmatrix}} - {{{\hat{h}}_{2}^{T}\left\lbrack {A_{r}^{T}\quad \gamma \quad S^{T}} \right\rbrack}\begin{pmatrix}Q \\{\gamma \quad d}\end{pmatrix}} - {{\left\lbrack {Q^{T}\quad \gamma \quad d^{T}} \right\rbrack \begin{bmatrix}A_{r} \\{\gamma \quad S}\end{bmatrix}}{\hat{h}}_{2}} +}} \\{\quad {{{{\hat{h}}_{2}^{T}\left\lbrack {A_{r}^{T}\quad \gamma \quad S^{T}} \right\rbrack}\begin{bmatrix}A_{r} \\{\gamma \quad S}\end{bmatrix}}{\hat{h}}_{2}}} \\{= \quad {{Q^{T}Q} + {\gamma^{2}d^{T}d} - {{\hat{h}}_{2}^{T}A_{r}^{T}Q} - {\gamma^{2}{\hat{h}}_{2}^{T}S^{T}d} - {Q^{T}A_{r}{\hat{h}}_{2}} - {\gamma^{2}d^{T}S{\hat{h}}_{2}} +}} \\{\quad {{{\hat{h}}_{2}^{T}A_{r}^{T}A_{r}{\hat{h}}_{2}} + {\gamma^{2}{\hat{h}}_{2}^{T}S^{T}F{\hat{h}}_{2}}}} \\{= \quad {{Q^{T}Q} - {{\hat{h}}_{2}^{T}A_{r}^{T}Q} - {q^{T}A_{r}{\hat{h}}_{2}} + {{\hat{h}}_{2}^{T}A_{r}^{T}A_{r}{\hat{h}}_{2}} +}} \\{\quad {\gamma^{2}\left\lbrack {{d^{T}d} - {{\hat{h}}_{2}^{T}F^{T}d} - {d^{T}{Fh}_{2}} + {{\hat{h}}_{2}^{T}F^{T}F{\hat{h}}_{2}}} \right\rbrack}} \\{= \quad {{{Q - {A_{r}{\hat{h}}_{2}}}}^{2} + {\gamma^{2}{{d - {S{\hat{h}}_{2}}}}^{2}}}}\end{matrix}$

[0058] where γ is the Lagrange multiplier and determines the relativeweighting the least squares solution of the reduced problem places inminimizing ∥Q−A_(r)ĥ₂∥² and ∥d−Fĥ₂∥². If γ is large, then the solutionwill minimize ∥d−Fĥ₂∥² more than ∥Q−A_(r)ĥ₂∥².

[0059] The smoothness of the least squares estimate ĥ₂ of the appendedreduced problem can be gauged by the norm of its second derivative withrespect to time according to: ${\hat{h}}_{2} = {{\begin{pmatrix}{{\hat{h}}_{2}(1)} \\{{\hat{h}}_{2}(2)} \\\vdots \\{{\hat{h}}_{2}\left( {N - M_{1}} \right)}\end{pmatrix}\frac{{\hat{h}}_{2}}{t}} = {{\frac{1}{\Delta \quad t}\begin{pmatrix}{{{\hat{h}}_{2}(1)} - {{\hat{h}}_{2}(0)}} \\{{{\hat{h}}_{2}(2)} - {{\hat{h}}_{2}(1)}} \\\vdots \\{{{\hat{h}}_{2}\left( {N - M_{1}} \right)} - {{\hat{h}}_{2}\left( {N - M_{1} - 1} \right)}}\end{pmatrix}} = {\frac{1}{\Delta \quad t}\begin{pmatrix}{{\hat{h}}_{2}(1)} \\{{\hat{h}(2)} - {\hat{h}(1)}} \\\vdots \\{{{\hat{h}}_{2}\left( {N - M_{1}} \right)} - {{\hat{h}}_{2}\left( {N - M_{1} - 1} \right)}}\end{pmatrix}}}}$$\frac{^{2}{\hat{h}}_{2}}{t^{2}} = {{\frac{1}{\Delta \quad t^{2}}\begin{pmatrix}{{{\hat{h}}_{2}(1)} - {2{{\hat{h}}_{2}(0)}} + {{\hat{h}}_{2}\left( {- 1} \right)}} \\{{{\hat{h}}_{2}(2)} - {2{{\hat{h}}_{2}(1)}} + {{\hat{h}}_{2}\left( {- 1} \right)}} \\\vdots \\{{{\hat{h}}_{2}\left( {N - M_{1}} \right)} - {2{{\hat{h}}_{2}\left( {N - M_{1} - 1} \right)}} + {{\hat{h}}_{2}\left( {N - M_{1} - 2} \right)}}\end{pmatrix}} = {{\frac{1}{\Delta \quad t^{2}}\begin{pmatrix}\begin{matrix}{{\hat{h}}_{2}(1)} \\{{{\hat{h}}_{2}(2)} - {2{{\hat{h}}_{2}(1)}}}\end{matrix} \\\vdots \\{{{\hat{h}}_{2}\left( {N - M_{1}} \right)} - {2{{\hat{h}}_{2}\left( {N - M_{1} - 1} \right)}} + {{\hat{h}}_{2}\left( {N - M_{1} - 2} \right)}}\end{pmatrix}} = {{\frac{1}{\Delta \quad t^{2}}\begin{bmatrix}1 & 0 & 0 & 0 & \cdots & 0 & 0 & 0 & 0 & 0 \\{- 2} & 1 & 0 & 0 & \cdots & 0 & 0 & 0 & 0 & 0 \\1 & {- 2} & 1 & 0 & \cdots & 0 & 0 & 0 & 0 & 0 \\0 & 1 & {- 2} & 1 & \cdots & 0 & 0 & 0 & 0 & 0 \\\quad & \quad & \quad & \quad & \cdots & \quad & \quad & \quad & \quad & \quad \\0 & 0 & 0 & 0 & \cdots & 0 & 1 & {- 2} & 1 & 0 \\0 & 0 & 0 & 0 & \cdots & 0 & 0 & 1 & {- 2} & 1\end{bmatrix}}\begin{pmatrix}{{\hat{h}}_{2}(1)} \\{{\hat{h}}_{2}(2)} \\{{\hat{h}}_{2}(3)} \\\vdots \\{{\hat{h}}_{2}\left( {N - M_{1} - 2} \right)} \\{{\hat{h}}_{2}\left( {N - M_{1} - 1} \right)} \\{\hat{h}\left( {N - M_{1}} \right)}\end{pmatrix}}}}$

[0060] To incorporate the smoothness constraint, the matrix S is setaccording to: $S = {\frac{1}{\Delta \quad t^{2}}\begin{bmatrix}1 & 0 & 0 & 0 & \cdots & 0 & 0 & 0 & 0 & 0 \\{- 2} & 1 & 0 & 0 & \cdots & 0 & 0 & 0 & 0 & 0 \\1 & {- 2} & 1 & 0 & \cdots & 0 & 0 & 0 & 0 & 0 \\0 & 1 & {- 2} & 1 & \cdots & 0 & 0 & 0 & 0 & 0 \\\quad & \quad & \quad & \quad & \cdots & \quad & \quad & \quad & \quad & \quad \\0 & 0 & 0 & 0 & \cdots & 0 & 1 & {- 2} & 1 & 0 \\0 & 0 & 0 & 0 & \cdots & 0 & 0 & 1 & {- 2} & 1\end{bmatrix}}$

[0061] where d=0. The least square solution of the appended reducedproblem facilitates minimizing both ∥Q−A_(r)ĥ₂∥² and ∥Sĥ₂∥² with therelative weighting controlled by γ.

[0062] In one embodiment, an inequality constraint is imposed on thesolution of the appended reduced problem. The appended reduced problemis: $\begin{pmatrix}Q \\{\gamma \quad d}\end{pmatrix} = {\begin{bmatrix}A_{r} \\{\gamma \quad S}\end{bmatrix}h_{2}}$

[0063] which can be written as: p=Eh₂.

[0064] In another embodiment, the least squares estimate of the appendedproblem under the linear inequality constraint Gh₂≧b is determined. Thedimension of the appended matrix E is M₂×(N−M₁). In one embodiment, Ehas a rank of K≦(N−M₁) and has the following singular valuedecomposition: $E = {{U\begin{bmatrix}S & 0 \\0 & 0\end{bmatrix}}V^{T}}$

[0065] where S is a K×K diagonal matrix with its diagonal entries equalto the singular values of E, U and V are orthogonal matrices ofdimensions M₂×M₂ and (N−M₁)×(N−M₁) respectively. The least squaressolution, ĥ₂, of the appended reduced problem p=Eh₂ will minimize∥Eh₂−p∥².

[0066] Using the singular value decomposition of E, ∥Eh₂−p∥² can besimplified as: $\begin{matrix}{{{{{Eh}_{2} - p}}^{2} = {{{{{U\begin{bmatrix}S & 0 \\0 & 0\end{bmatrix}}V^{T}h_{2}} - p}}^{2} = {{{U^{T}{U\begin{bmatrix}S & 0 \\0 & 0\end{bmatrix}}V^{T}h_{2}} - {U^{T}p}}}^{2}}},{{Since}\quad U^{T}\quad {is}\quad {orthogonal}}} \\{= {{{{\begin{bmatrix}S & 0 \\0 & 0\end{bmatrix}V^{T}h} - {U^{T}p}}}^{2} = {{{{{\begin{bmatrix}S & 0 \\0 & 0\end{bmatrix}\begin{bmatrix}V_{1}^{T} \\V_{2}^{T}\end{bmatrix}}h_{2}} - {\begin{bmatrix}U_{1}^{T} \\U_{2}^{T}\end{bmatrix}p}}}^{2} = {{\begin{bmatrix}{{SV}_{1}^{T}h_{2}} \\0\end{bmatrix} - \begin{bmatrix}{U_{1}^{T}p} \\{U_{2}^{T}p}\end{bmatrix}}}^{2}}}} \\{= {{\begin{bmatrix}{{{SV}_{1}^{T}h_{2}} - {U_{1}^{T}p}} \\{{- U_{2}^{T}}p}\end{bmatrix}}^{2} = {{{{{SV}_{1}^{T}h_{2}} - {U_{1}^{T}p}}}^{2} + {{U_{2}^{T}p}}^{2}}}}\end{matrix}$

[0067] Since ∥U₂ ^(T)p∥² is constant, minimizing ∥Eh₂−p∥² is equivalentto minimizing ∥SV₁ ^(T)h₂−U₁ ^(T)P∥².

[0068] In one embodiment, the variable h₂ is changed to z such thatz=SV₁ ^(T)h₂−U₁ ^(T)p

V₁ ^(T)h₂=S⁻¹z+S⁻¹U₁ ^(T)p∥².

[0069] Minimizing ∥Eh₂−p∥² is the same as minimizing ∥z∥² and, thereforethe inequality constraint becomes:

Gh₂≧b

G(V ₁ V ₁ ^(T) +V ₂ V ₂ ^(T))h ₂ ≧b

GV ₁(S ⁻¹ z+S ⁻¹ U ₁ ^(T) p)+GV ₂ V ₂ ^(T) h ₂ ≧b

GV ₁(S ⁻¹ z+S ⁻¹ U ₁ ^(T) p)+GV ₂ V ₁ ^(T) V ₁(S ⁻¹ z+S ⁻¹ U ₁ ^(T) p)≧b

GV ₁ S ⁻¹ z+GV ₂ V ₂ ^(T) S ⁻¹ z≧b−GV ₁ S ⁻¹ U ₁ ^(T) p−GV ₂ V ₂ ^(T) V₁ S ⁻¹ U ₁ ^(T) p

G(I ₂ +V ₂ V ₂ ^(T))V ₁ S ⁻¹ z≧b−G(I ₂ +V ₂ V ₂)V ₁ S ⁻¹ U ₁ ^(T) p

[0070] Therefore, the problem of minimizing ∥Eh₂=p∥² subject to theconstraint Gh₂≧b is equivalent to minimizing ∥z∥² subject to G(I₂+V₂V₂^(T))V₁S⁻¹z≧b−G(I₂+V₂V₂ ^(T))V₁S⁻¹U₁ ^(T)p. From the solution,{circumflex over (z)}, that minimize ∥z∥², the corresponding leastsquares solution, ĥ₂, of the appended reduced problem Eh₂=p can be foundby the inversion of the equation z=SV₁ ^(T)h₂−U₁ ^(T)p. From the leastsquares solution, ĥ₂, of the appended reduced problem, Eh₂=p, the fullimpulse residue function, h, can be reconstituted by:$h = \begin{pmatrix}h_{1} \\{\hat{h}}_{2}\end{pmatrix}$

[0071] and h₁ is a (I+L−1)×1 vector whose first I elements are zero andthe subsequent L−1 elements are equal to the first element of ĥ₂.

[0072] In another embodiment, E has a rank of K=(N−M₁) and has thefollowing singular value decomposition: ${E = {{U\begin{bmatrix}S \\0\end{bmatrix}}V^{T}}},$

[0073] where S is a K×K diagonal matrix with its diagonal entries equalto the singular values of E, U and V are orthogonal matrices ofdimensions M₂×M₂ and (N−M₁)×(N−M₁) respectively. The least squaressolution, ĥ₂, of the appended reduced problem f=Eh₂ facilitatesminimizing ∥Eh₂−p∥².

[0074] Using a singular value decomposition of E, ∥Eh₂−p∥² can besimplified as: $\begin{matrix}{{{{{Eh}_{2} - p}}^{2} = {{{{{U\begin{bmatrix}S \\0\end{bmatrix}}V^{T}h_{2}} - p}}^{2} = {{{U^{T}{U\begin{bmatrix}S \\0\end{bmatrix}}V^{T}h_{2}} - {U^{T}p}}}^{2}}},{{Since}\quad U^{T}\quad {is}\quad {orthogonal}}} \\{= {{{{\begin{bmatrix}S \\0\end{bmatrix}V^{T}h_{2}} - {U^{T}p}}}^{2} = {{{{\begin{bmatrix}S \\0\end{bmatrix}V^{T}h_{2}} - {\begin{bmatrix}U_{1}^{T} \\U_{2}^{T}\end{bmatrix}p}}}^{2} = {{\begin{bmatrix}{{SV}^{T}h_{2}} \\0\end{bmatrix} - \begin{bmatrix}{U_{1}^{T}p} \\{U_{2}^{T}p}\end{bmatrix}}}^{2}}}} \\{= {{\begin{bmatrix}{{{SV}^{T}h_{2}} - {U_{1}^{T}p}} \\{{- U_{2}^{T}}p}\end{bmatrix}}^{2} = {{{{{SV}^{T}h_{2}} - {U_{1}^{T}p}}}^{2} + {{U_{2}^{T}p}}^{2}}}}\end{matrix}$

[0075] since ∥U₂ ^(T)p∥² is constant, minimizing ∥Eh₂−p∥² is equivalentto minimizing ∥SV^(T)h₂−U₁ ^(T)p∥². In one embodiment, the variable h₂is changed to z such that: z=SV^(T)h₂−U₁ ^(T)p

V^(T)h₂=S⁻¹z+S⁻¹U₁ ^(T)p

[0076] Minimizing ∥Eh₂−p∥² is the same as minimizing ∥z∥² and theinequality constraint becomes:

Gd₂≧b

G(VV ₁ ^(T))h ₂ ≧b

GV(S ⁻¹ z+S ⁻¹ U ₁ ^(T) p)≧b

GVS ⁻¹ z≧b−GVS ⁻¹ U ₁ ^(T) p

[0077] Minimizing ∥Eh₂−p∥² subject to the constraint Gh₂≧b is equivalentto minimizing ∥z∥² subject to GVS⁻¹z≧b−GVS⁻¹U₁ ^(T)p.

[0078] A formulation of smoothness, nonnegative and monotonicityconstraints for the reduced problem can now be stated. The least squaressolution of the reduced problem: Q=A_(r)h₂ with the smoothness,nonnegative and monotonicity constraints can be recast as the leastsquares solution of the appended reduced problem according to:$\begin{pmatrix}q \\{\gamma \quad d}\end{pmatrix} = {{\begin{bmatrix}A_{r} \\{\gamma \quad S}\end{bmatrix}h_{2}\quad {or}\quad p} = {Eh}_{2}}$

[0079] with the constraint I_((N−M) ₁ _()×(N−M) ₁ ₎h₂≧0 and, Dh₂≧0.where D is a (N−M₁−1)×(N−M₁) matrix of the form: $D = \begin{bmatrix}1 & {- 1} & 0 & \cdots & \cdots & \cdots & 0 & 0 & 0 \\0 & 1 & {- 1} & \cdots & \cdots & \cdots & 0 & 0 & 0 \\\vdots & \vdots & \vdots & \vdots & \vdots & \vdots & \vdots & \vdots & \vdots \\\vdots & \vdots & \vdots & \vdots & \vdots & \vdots & \vdots & \vdots & \vdots \\0 & 0 & 0 & \cdots & \cdots & \cdots & 1 & {- 1} & 0 \\0 & 0 & 0 & \cdots & \cdots & \cdots & 0 & 1 & {- 1}\end{bmatrix}$

[0080] such that Dh₂≧0 is equivalent to h₂(1)≧h₂(2)≧ . . . h₂(N−M₁).

[0081] The two inequality constraints I_((N−M) ₁ _()×(N−M) ₁ ₎h₂≧0 andDh₂≧0 can be combined into ${\begin{bmatrix}I_{{({N - M_{1}})}{x{({N - M_{1}})}}} \\D_{{({N - M_{1} - 1})}{x{({N - M_{1}})}}}\end{bmatrix}h_{2}} \geq 0$

[0082] or, Gh₂≧0 where G is a 2(N−M₁1)×(N−M₁) matrix as defined hereinsuch that I_((N−M) ₁ _()×(N−M) ₁ ₎ is an (N−M₁)×(N−M₁) identity matrix.

[0083] The least squares solution of the appended reduced problem p=Eh₂subject to the nonnegative and montonicity constraints can be determinedusing the method described herein. From the least squares solution, ĥ₂,of the appended reduced problem, Eh₂=p, the full impulse residuefunction, h, can be reconstituted by: $h = \begin{pmatrix}h_{1} \\{\hat{h}}_{2}\end{pmatrix}$

[0084] and h₁ is a (I+L−1)×1 vector whose first I elements are zero andthe A subsequent L−1 elements are equal to the first element of ĥ₂.

[0085] In one embodiment, in tissue where there is no leakage ofcontrast from the blood stream into the interstitial space, for example,the brain with intact blood-brain-barrier, the reconstituted fullimpulse residue function h(t) can be used to determine Tissue Blood Flow(TBF) as the peak height of h(t), Tissue Blood Volume (TBV) as the areaunderneath h(t) and Tissue Mean Transit Time (TMTT) as the TBV dividedby TBF.

[0086] In one embodiment, in tissue where there is leakage of contrastfrom the blood stream into the interstitial space, the impulse residuefunction h(t) can be modeled by the Johnson-Wilson model.

[0087] An adiabatic approximation of the Johnson-Wilson Model is thenused to parameterize the impulse residue function in terms of TBF, TBV,TMTT and TPS. In one embodiment, the impulse residue function for theadiabatic approximation of the Johnson-Wilson model is: $\begin{matrix}{{R(t)} = \left\{ \begin{matrix}0.0 & {t \leq t_{o}} \\1.0 & {t_{o} \leq t \leq {t_{o} + W}} \\{Ee}^{- {k{({t - t_{o} - W})}}} & {t \geq {t_{o} + W}}\end{matrix} \right.} & (1)\end{matrix}$

[0088] where t_(o) is a time delay between an arterial curve of contrastconcentration, C_(a)(t), and the tissue residue function, Q(t), W is thetissue mean transit time (TMTT), and k is the rate constant of theexponential function and is defined according to: $\begin{matrix}{k = \frac{FE}{V_{e}}} & (2)\end{matrix}$

[0089] where F is the Tissue Blood Flow (TBF), E is the extractionefficiency of contrast material from blood and V_(e) is the distributionspace of contrast material in the intersitial space. Tissue Blood Volumeis equal to TBF×TMTT. Tissue capillary permeability surface area product(TPS) is: −ln(1−E).

[0090] The tissue residue function, Q(t), is given by:

Q(t)=F·[C _(a)(t)*R(t)]  (3)

[0091] where * is the convolution operation.

[0092] The linearization of the tissue residue function (or theexpression of the tissue residue function or some t function of it interms of linear functions of the parameters (or some combination ofthem) of the model) is based on finding the (time) integral of Q(t)according to: $\begin{matrix}{{\int_{0}^{T}{{Q(t)}{t}}} = {\int_{0}^{T}\quad {{t}{\int_{0}^{t}{{C_{a}(u)}\quad {R\left( {t - u} \right)}{u}}}}}} & \left( {4a} \right)\end{matrix}$

[0093] Interchanging the order of integration: $\begin{matrix}{{\int_{0}^{T}{{Q(t)}{t}}} = {{\int_{0}^{T}\quad {{u}{\int_{u}^{T}{{C_{a}(u)}\quad {R\left( {t - u} \right)}{t}}}}} = {{\int_{0}^{T}{{C_{a}(u)}{u}{\int_{u}^{T}{{R\left( {t - u} \right)}{t}}}}} = {\int_{0}^{T}{{C_{a}(u)}{u}{\int_{0}^{T - u}{{R(t)}{t}}}}}}}} & \left( {4b} \right)\end{matrix}$

[0094] The time integral ∫₀^(T)Q(t)t

[0095] is evaluated for the three time intervals

[0096] (iii) (i) 0≦T≦t_(o)

[0097] (iv) (ii) t_(o)≦T≦t_(o)+W

[0098] (V) t_(o)+W≦T

[0099] For example when the time integral ∫₀^(T)Q(t)t

[0100] is evaluated for the time interval 0≦T≦t_(o), then∫₀^(T)Q(t)t = ∫₀^(T)C_(a)(u)u∫₀^(T − u)R(t)t

[0101] for 0≦u≦T

0≧−u≧−T

T≧T−u≧0

0≦T−u≦T≦t_(o), therefore: $\begin{matrix}{{\int_{0}^{T}{{Q(t)}{t}}} = 0} & (5)\end{matrix}$

[0102] When the time integral ∫₀^(T)Q(t)t

[0103] is evaluated for the time interval t_(o)≦T≦t_(o)+W, then settingT=t_(o)+a, where a≦W, the time integral is:∫₀^(T)Q(t)t = F∫₀^(t_(o) + a)C_(a)(u)u∫₀^(t_(o) + a − u)R(t)t

$= {{{F \cdot \left\lbrack {{\int_{0}^{a}{{C_{a}(u)}{u}}} + {\int_{a}^{t_{o} + a}{{C_{a}(u)}{u}}}} \right\rbrack}\left( {\int_{0}^{t_{o} + a - u}{{R(t)}{t}}} \right)} = {{F \cdot {\int_{o}^{a}{{C_{a}(u)}{u_{\underset{A}{}}}{\int_{0}^{t_{o} + a - u}{{R(t)}{t}}}}}} + {F \cdot {\int_{a}^{t_{o} + a}{{C_{a}(u)}{u}{\int_{0}^{t_{o} + a - u}{{R(t)}{t}}}}}}}}$

[0104] A: For 0≦u≦a,t_(o)+a≦t_(o)+a−u≧t_(o) or t_(o)+W≧t_(o)+a−u≧t_(o),thenF∫₀^(a)C_(a)(u)u∫₀^(t_(o) + a − u)R(t)t = F∫₀^(a)C_(a)(u)u[∫_(o)^(t_(o))R(t)t + ∫_(t_(o))^(t_(o) + a − u)R(t)t] = F_(′)∫₀^(a)C_(a)(u)u∫_(t_(o))^(t_(o) + a − u)R(t)t = F(a − u)∫₀^(a)C_(a)(u)u.

[0105] B: For a≦u≦a+t_(o),t_(o)≧t_(o)+a−u≧0, thenF * ∫_(a)^(t_(o) + a)C_(a)(u)  u∫₀^(t_(o) + a − u)R(t)  t = 0

[0106] Combining A and B above: $\begin{matrix}{{{\int_{0}^{T}{{Q(t)}\quad {t}}} = {F*{\int_{0}^{a}{\left( {a - u} \right){C_{a}(u)}\quad {u}\quad {or}}}}},\quad {{F\left\lbrack {{\left( {T - t_{o}} \right){\int_{0}^{T - t_{o}}{{C_{a}(u)}\quad {u}}}} - {\int_{0}^{T - t_{o}}{{{uC}_{a}(u)}\quad {u}}}} \right\rbrack} = {\int_{0}^{T}{{Q(t)}\quad {t}}}}} & \left( {6a} \right)\end{matrix}$

[0107] When the time integral ∫₀^(T)Q(t)  t

[0108] is evaluated for the time interval T≧t_(o)+W, then and settingT=t_(o)+W+a, where a≧0 then${\int_{0}^{T}{{Q(t)}\quad {t}}} = {{F \cdot {\int_{0}^{T}{{C_{a}(u)}\quad {u}{\int_{0}^{T - u}{{R(t)}\quad {t}}}}}}\quad  = {{F*{\int_{0}^{t_{o} + a + W}{{C_{a}(u)}\quad {u}{\int_{0}^{t_{o} + a + W - u}{{R(t)}\quad {t}}}}}}\quad  = {{F\left\lbrack {{\int_{0}^{a}{{C_{a}(u)}\quad {u}}} + {\int_{a}^{a + W}{{C_{a}(u)}\quad {u}}} + {\int_{a + W}^{t_{o} + a + W}{{C_{a}(u)}\quad {u}{\int_{0}^{t_{o} + a + W - u}{{R(t)}\quad {t}}}}}} \right\rbrack} = {\underset{\underset{C}{}}{F{\int_{0}^{a}{{C_{a}(u)}\quad {u}{\int_{0}^{{t\quad 0} + a + w - u}{{R(t)}\quad {t}}}}}} + \underset{\underset{D}{}}{F{\int_{0}^{a + W}{{C_{a}(u)}\quad {u}{\int_{0}^{{t\quad 0} + a + W - u}{{R(t)}\quad {t}}}}}} + \underset{\underset{E}{}}{F{\int_{a + W}^{{t\quad 0} + a + W}{{C_{a}(u)}\quad {u}{\int_{0}^{{t\quad 0} + a + W - u}{{R(t)}\quad {t}}}}}}}}}}$

[0109] C: For 0≦u≦a,t_(o)+a+W≧t_(o)+a+W−u≧t_(o)+W, then${{{{{{{{F{\int_{0}^{a}{{C_{a}(u)}\quad {u}{\int_{0}^{t_{o} + a + W - u}{{R(t)}\quad {t}}}}}} = {F{\int_{0}^{a}{{C_{a}(u)}\quad {u}}}}}\quad}{\quad{\left\lbrack {{\int_{0}^{t_{o}}{{R(t)}\quad {t}}} + {\int_{0}^{t_{o} + W}{{R(t)}\quad {t}}} + {\int_{t_{o} + W}^{t_{o} + a + W - u}{{R(t)}\quad {t}}}} \right\rbrack  = {{F{\int_{0}^{a}{{C_{a}(u)}\quad {{u\left\lbrack {W + {\int_{t_{o} + W}^{t_{o} + a + W - u}{{R(t)}\quad {t}}}} \right\rbrack}}}}} = {{{{FW}{\int_{0}^{a}{{C_{a}(u)}\quad {u}}}} + {F{\int_{0}^{a}{{C_{a}(u)}\quad {u}{\int_{t_{o} + W}^{t_{o} + a + W - u}{{R(t)}\quad {t}}}}}}} = {{{{FW}{\int_{0}^{a}{{C_{a}(u)}\quad {u}}}} + {F{\int_{0}^{a}{{C(u)}\quad {u}{\int_{t_{o} + W}^{t_{o} + a + W - u}{E\quad ^{- {k{({t - t_{o} - W})}}}\quad {t}}}}}}} = {{{{FW}{\int_{0}^{a}{{C_{a}(u)}\quad {u}}}} + {{FE}{\int_{0}^{a}{{C_{a}(u)}\quad {u}{\int_{0}^{a - u}{^{- {kt}}\quad {t}}}}}}} = {{{{FW}{\int_{0}^{a}{{C_{a}(u)}\quad {u}}}} - {\frac{{FE}^{a}}{k}{\int_{0}^{a}{{C_{a}(u)}\quad {{u\left\lbrack {^{- {k{({a - u})}}} - 1} \right\rbrack}}}}}}\quad  = {{{FW}{\int_{0}^{a}{{C_{a}(u)}\quad {u}}}} +}}}}}}}\quad}{\quad{{{{{\frac{{FE}^{a}}{k}{\int_{0}^{a}{{C_{a}(u)}\quad {u}}}} - \underset{\underset{F}{}}{\frac{{FE}^{a}}{k}{\int_{0}^{a}{{C_{a}(u)}\quad ^{- {k{({a - u})}}}{u}}}}}\quad {Q(T)}} = {{{FW}{\int_{0}^{T}{{C_{a}(u)}\quad {R\left( {T - u} \right)}{u}}}} = {{{F{\int_{0}^{t_{o} + a + W}{{C_{a}(u)}{R\left( {t_{o} + a + W - u} \right)}\quad {u}}}}{F:}}\quad  = {{{F{\int_{0}^{a}{{C_{a}(u)}\quad {R\left( {t_{o} + a + W - u} \right)}{u}}}} + {F{\int_{a}^{a + W}{{C_{a}(u)}{R\left( {t_{o} + a + W - u} \right)}\quad {u}}}} + {F{\int_{a + W}^{t_{o} + a + W}{{C_{a}(u)}{R\left( {t_{o} + a + W - u} \right)}\quad {u}}}}} = {{{FE}{\int_{0}^{a}{{C_{a}(u)}^{- {k{({a - u})}}}\quad {u}}}} + {F{\int_{a}^{a + W}{{C_{a}(u)}\quad {u}\quad {Therefore}}}}}}}}},}\quad\quad}\frac{FE}{k}{\int_{0}^{a}{{C_{a}(u)}^{- {k{({a - u})}}}\quad {u}}}} =}\quad}\frac{Q(T)}{k}} - {\frac{F}{k}{\int_{a}^{a + W}{{C_{a}(u)}\quad {u}}}}}\quad$and  ${F{\int_{0}^{a}{{C_{a}(u)}\quad {u}{\int_{0}^{t_{o} + a + W - u}{{R(t)}\quad {t}}}}}} = {{{FW}{\int_{0}^{a}{{C_{a}(u)}\quad {u}}}} + {\frac{FE}{k}{\int_{0}^{a}{{C_{a}(u)}\quad {u}}}} - \frac{Q(T)}{k} + {\frac{F}{k}{\int_{a}^{a + W}{{C_{a}(u)}\quad {u}}}}}$C:  

[0110] D: For a≦u≦a+W,t_(o)+W−u≧t_(o)+a+W−u≧t_(o)

t_(o)+W≧t_(o)+a+W−u≧t_(o)F∫_(a)^(a + W)C_(a)(u)  u∫₀^(t_(o) + a + W − u)R(t)  t = F∫_(a)^(a + W)C_(a)(u)  u[∫_(a)^(t_(o))R(t)  t + ∫_(t_(o))^(t_(o) + a + W − u)R(t)  t] = F∫_(a)^(a + W)C_(a)(u)  u∫_(t_(o))^(t_(o) + a + W − u)R(t)  t = F[a + W − u]∫_(a)^(a + W)C_(a)(u)  u  

[0111] E: For a+W≦u≦a+W+t_(o),t_(o)≧t_(o)+a+W−u≧0F∫_(a + W)^(t_(o) + a + W)C_(a)(u)  u∫₀^(t_(o) + a + W − u)R(t)  t = 0

[0112] Therefore, for T−t_(o)+W+a, where a≧0: $\begin{matrix}{{{\int_{0}^{T}{{Q(t)}\quad {t}}} = {{{FW}{\int_{0}^{a}{{C_{a}(u)}\quad {u}}}} + {\frac{FE}{k}{\int_{0}^{a}{{C_{a}(u)}\quad {u}}}} - \frac{Q(T)}{k} + {\frac{F}{k}{\int_{a}^{a + W}{{C_{a}(u)}\quad {u}}}} + {{F\left\lbrack {a + W - u} \right\rbrack}{\int_{0}^{a + W}{{C_{a}(u)}\quad {u}}}}}}{{or},\quad {{QT} = {{{Fk}\left\lbrack {{W{\int_{0}^{a + W}{{C_{a}(u)}\quad {u}}}} + {a{\int_{a}^{a + W}{{C_{a}(u)}\quad {u}}}} - {\int_{a}^{a + W}{{{uC}_{a}(u)}\quad {u}}}} \right\rbrack} + {F{\int_{a}^{a + W}{{C_{a}(u)}\quad {u}}}} + {{FE}{\int_{0}^{a + W}{{C_{a}(u)}\quad {u}}}} - {k{\int_{0}^{T}{{Q(t)}\quad {t}}}}}}}} & \left( {7a} \right) \\{= {{{Fk}\left\lbrack {{W{\int_{0}^{T + t_{o}}{{C_{a}(u)}\quad {u}}}} + {\left( {T - t_{o} - W} \right){\int_{T - t_{o} - W}^{T - t_{o}}{{C_{a}(u)}\quad {u}}}} - {\int_{T - t_{o} - W}^{T - t_{o}}{{{uC}_{a}(u)}\quad {u}}}} \right\rbrack} + {F{\int_{T - t_{o} - W}^{T - t_{o}}{{C_{a}(u)}\quad {u}}}} + {{FE}{\int_{0}^{T - t_{o} - W}{{C_{a}(u)}\quad {u}}}} - {k{\int_{0}^{T}{{Q(t)}\quad {t}}}}}} & \left( {7b} \right)\end{matrix}$

[0113] In one embodiment, the time integral ∫₀^(T)Q(t)  t

[0114] for the three time intervals

[0115] (i) 0≦T≦t_(o)

[0116] (ii) t_(o)≦T≦t_(o)+W

[0117] (iii) t_(o)+W≦T

[0118] generates the following system of linear equations0 = −k∫₀^(T)Q(t)  t  0 ≤ T ≤ t_(o)0 = Fk[(T − t_(o))∫₀^(T − t_(o))C_(a)(u)  u − ∫₀^(T − t_(o))uC_(a)(u)  u] − k∫₀^(T)Q(t)  t  t_(o) ≤ T ≤ t_(o) + WQ(T) = Fk[W∫₀^(T − t_(o))C_(a)(u)  u + (T − t_(o) − W)∫_(T − t_(o) − W)^(T − t_(o))C_(a)(u)  u − ∫_(T − t_(o) − W)^(T − t_(o))uC_(a)(u)  u] + F∫_(T − t_(o) − W)^(T − t_(o))C_(a)(u)  u + FE∫₀^(T − t_(o) − W)C_(a)(u)  u −     k∫₀^(T)Q(t)  tT   ≥   t_(o) +   W

[0119] involving a plurality of unknowns Fk, F, FE and k, provided thatt₀ and W are known.

[0120] For a set of given t₀ and W, a least squares algorithm is used toestimate Fk, F, FE and k subject to the following linear constraints:

[0121] Fk≧0

[0122] F≧0

[0123] FE≧0

[0124] k≧1

[0125] Fk≧F

[0126] F≧FE

[0127] The algorithm for estimating t₀, W, Fk, F, FE and k includes twomain golden section search subroutines. An outer golden section searchsubroutine searches for t₀, an inner golden section search subroutinesearches for W with the value of t₀, assumed in the outer routine.Within the inner golden section search routine, values of t₀ and W arefixed and the optimization of Fk, F, FE and k can proceed with LDPalgorithms.

[0128] In one embodiment, values of F, E, Ve and k are: F=0.1 ml/min/g,E=0.5, Ve0.25, interstitial space, and k=FE/Ve=0.2 min−1=12s−1. Inanother embodiment, values of F, E, V_(e) and k are: F=0.1 ml/min/g,E=0.25, V_(e)0.25, interstitial space, and k=FE/V_(e)=0.1 min⁻¹=6s⁻¹.

[0129] Although k>=1 is the constraint, after solution, for any k>5s⁻¹assumes that there is no leakage and do not calculate PS, andrecalculate blood volume and MTT instead of taking W as MTT and EW asvolume. The volume can be taken as the area underneath the flow scaledimpulse residue function (F.R(t)), and MTT as the volume divided by theheight of the flow scaled impulse residue function.

[0130] Discretization of the system of linear equations is performedaccording to:

[0131] Let Δt be the time interval, and t₀=mΔt, W=nΔt, and T=iΔt, then:$\begin{matrix}\begin{matrix}{0 = \quad {{{- k}{\sum\limits_{j = 1}^{i}\quad {{Q(j)}\Delta \quad t\quad 0}}} \leq i \leq m}} \\{0 = \quad {{{{Fk}\left\lbrack {{\left( {i - m} \right)\Delta \quad t{\sum\limits_{j = 1}^{i - m}\quad {{C_{a}(j)}\Delta \quad t}}} - {\sum\limits_{j = 1}^{i - m}\quad {j\quad {C_{a}(j)}\Delta \quad t}}} \right\rbrack} - {k{\sum\limits_{j = 1}^{i}\quad {{Q(j)}\Delta \quad {tm}}}}} \leq i \leq {n + m}}} \\{{Q(T)} = \quad {{{Fk}\left\lbrack {{W{\sum\limits_{j = 1}^{i - m}\quad {{C_{a}(j)}\Delta \quad t}}} + {\left( {i - m - n} \right)\Delta \quad t{\sum\limits_{j = {i - m - n}}^{i - m}\quad {{C_{a}(j)}\Delta \quad t}}} - {\sum\limits_{j = {i - m - n}}^{i - m}\quad {{C_{a}(j)}\Delta \quad t}}} \right\rbrack} +}} \\{\quad {{{F{\sum\limits_{j = {i - m - n}}^{i - m - n}\quad {{C_{a}(j)}\Delta \quad t}}} + {{FE}{\sum\limits_{j = 1}^{i - m - n}\quad {{C_{a}(j)}\Delta \quad t}}} - {k{\sum\limits_{j = 1}^{i}\quad {{Q(j)}\Delta \quad {ti}}}}} \geq {n + m}}}\end{matrix} & (9)\end{matrix}$

[0132] The various sums involving the arterial function C_(a)(t) can bepre-calculated and stored in 2-dimensional arrays AAR and AFM. Forexample,$\sum\limits_{j = {i - m - n}}^{i - m}\quad {{C_{a}(j)}\Delta \quad t}$

[0133] as the [(i−m−n),(i−m)]th element of Area array referred to hereinas an arterial area array (AAR). Also,$\sum\limits_{j = {i - m - n}}^{i - m}\quad {{C_{a}(j)}\Delta \quad t}$

[0134] as the [(i−m−n),(i−m)]th element of First-moment array referredto herein as an arterial first moment array (AFM). The various sumsinvolving the tissue residue function Q(t) can be pre-calculated andstored in a 1-dimensional array AQ. For example,$\sum\limits_{j = 1}^{i - m}\quad {{Q(j)}\Delta \quad t}$

[0135] as the ith element of an array for area of tissue residuefunction Q(t), (AQ).

[0136] A solution for F, E, and k when to and W are known can bedetermined as follows:

[0137] Using the arrays: Q, AAR, AFM and AQ, Eq. (9) can be written asEq. (10): $\begin{bmatrix}{{Q(1)} = 0} \\\vdots \\{{Q(m)} = 0} \\{{Q\left( {m + 1} \right)} = 0} \\{{Q\left( {m + 2} \right)} = 0} \\\vdots \\{Q\left( {M + n} \right)} \\{Q\left( {m + n + 1} \right)} \\{Q\left( {m + n + 2} \right)} \\\vdots \\{Q\left( {m + n + p} \right)}\end{bmatrix} = {\begin{bmatrix}{\quad 0} & 0 & 0 & {- {{AQ}(1)}} \\{\quad \vdots} & \vdots & \vdots & \vdots \\{\quad 0} & 0 & 0 & {- {{AQ}(m)}} \\{\quad {{\Delta \quad {t \cdot {{AAR}\left( {1,1} \right)}}} - {{AFM}\left( {1,1} \right)}}} & 0 & 0 & {- {{AQ}\left( {m + 1} \right)}} \\{\quad {{2\Delta \quad {t \cdot {{AAR}\left( {1,2} \right)}}} - {{AFM}\left( {1,2} \right)}}} & 0 & 0 & {- {{AQ}\left( {m + 2} \right)}} \\\vdots & \vdots & \vdots & \vdots \\{\quad {{{n \cdot \Delta}\quad {t \cdot {{AAR}\left( {1,n} \right)}}} - {{AFM}\left( {1,n} \right)}}} & 0 & 0 & {- {{AQ}\left( {m + n} \right)}} \\{\quad {{n \cdot \Delta \cdot {{AAR}\left( {1,{n + 1}} \right)}} + {\Delta \quad {t \cdot {{AAR}\left( {1,{n + 1}} \right)}}} - {{AFM}\left( {1,{n + 1}} \right)}}} & {{AAR}\left( {1,1} \right)} & {{AAR}\left( {1,1} \right)} & {- {{AQ}\left( {m + n + 1} \right)}} \\{{{n \cdot \Delta}\quad {t \cdot {{AAR}\left( {1,{n + 2}} \right)}}} + {{2 \cdot \Delta}\quad {T \cdot {{AAR}\left( {2,{n + 2}} \right)}}} - {{AFM}\left( {2,{n + 2}} \right)}} & {{AAR}\left( {2,{n + 2}} \right)} & {{AAR}\left( {1,2} \right)} & {- {{AQ}\left( {m + n + 2} \right)}} \\\vdots & \vdots & \vdots & \vdots \\{{{n \cdot \Delta}\quad {t \cdot {{AAR}\left( {1,{n + p}} \right)}}} + {{p \cdot \Delta}\quad {t \cdot {{AAR}\left( {p,{n + p}} \right)}}} - {{AFM}\left( {p,{n + p}} \right)}} & {{AAR}\left( {p,{n + p}} \right)} & {{AAR}\left( {1,p} \right)} & {- {{AQ}\left( {m + n + p} \right)}}\end{bmatrix}{\quad\begin{bmatrix}{Fk} \\F \\{FE} \\k\end{bmatrix}}}$

[0138] or Q=Ax

[0139] and: $Q = {{\begin{bmatrix}{{Q(1)} = 0} \\\vdots \\{{Q(m)} = 0} \\{{Q\left( {m + 1} \right)} = 0} \\{{Q\left( {m + 2} \right)} = 0} \\\vdots \\{Q\left( {M + n} \right)} \\{Q\left( {m + n + 1} \right)} \\{Q\left( {m + n + 2} \right)} \\\vdots \\{Q\left( {m + n + p} \right)}\end{bmatrix}A} = {{\begin{bmatrix}{\quad 0} & 0 & 0 & {- {{AQ}(1)}} \\{\quad \vdots} & \vdots & \vdots & \vdots \\{\quad 0} & 0 & 0 & {- {{AQ}(m)}} \\{\quad {{\Delta \quad {t \cdot {{AAR}\left( {1,1} \right)}}} - {{AFM}\left( {1,1} \right)}}} & 0 & 0 & {- {{AQ}\left( {m + 1} \right)}} \\{\quad {{2\Delta \quad {t \cdot {{AAR}\left( {1,2} \right)}}} - {{AFM}\left( {1,2} \right)}}} & 0 & 0 & {- {{AQ}\left( {m + 2} \right)}} \\\vdots & \vdots & \vdots & \vdots \\{\left. \quad {{{n \cdot \Delta \cdot \Delta}\quad {t \cdot {AA}}\quad 1},n} \right) - {{AFM}\left( {1,n} \right)}} & 0 & 0 & {- {{AQ}\left( {m + n} \right)}} \\{\left. \quad {{{n \cdot \Delta \cdot \Delta \cdot {AA}}\quad 1},{n + 1}} \right) + {\Delta \quad {t \cdot {{AAR}\left( {1,{n + 1}} \right)}}} - {{AFM}\left( {1,{n + 1}} \right)}} & {{AAR}\left( {1,{n + 1}} \right)} & {{AAR}\left( {1,1} \right)} & {- {{AQ}\left( {m + n + 1} \right)}} \\{\left. {{\left. {{{n \cdot \Delta \cdot \Delta}\quad {t \cdot {AA}}\quad 1},{n + 2}} \right) + {{2 \cdot {\Delta\Delta}}\quad {T \cdot {AAR}}\quad 2}},{n + 2}} \right) - {{AFM}\left( {2,{n + 2}} \right)}} & {{AAR}\left( {2,{n + 2}} \right)} & {{AAR}\left( {1,2} \right)} & {- {{AQ}\left( {m + n + 2} \right)}} \\\vdots & \vdots & \vdots & \vdots \\{\left. {\left. {{{n \cdot \Delta \cdot \Delta}\quad {t \cdot {AA}}\quad 1},{n + p}} \right) + {{p \cdot \Delta \cdot \Delta}\quad {t \cdot {AARn}}} + p} \right) - {{AFM}\left( {p,{n + p}} \right)}} & {{AAR}\left( {p,{n + p}} \right)} & {{AAR}\left( {1,p} \right)} & {- {{AQ}\left( {m + n + p} \right)}}\end{bmatrix}x} = {\quad\begin{bmatrix}{Fk} \\F \\{FE} \\k\end{bmatrix}}}}$

[0140] The constraints are: $\begin{matrix}{{{{\begin{bmatrix}1 & 0 & 0 & 0 \\0 & 1 & 0 & 0 \\0 & 0 & 1 & 0 \\0 & 0 & 0 & 1 \\1 & {- 1} & 0 & 0 \\0 & 1 & {- 1} & 0\end{bmatrix}\begin{bmatrix}{Fk} \\F \\{FE} \\k\end{bmatrix}} = {\begin{bmatrix}0 \\0 \\0 \\1 \\0 \\0\end{bmatrix}\quad {or}}},{{Gx} \geq b}}{G = {{\begin{bmatrix}1 & 0 & 0 & 0 \\0 & 1 & 0 & 0 \\0 & 0 & 1 & 0 \\0 & 0 & 0 & 1 \\1 & {- 1} & 0 & 0 \\0 & 1 & {- 1} & 0\end{bmatrix}\quad {and}\quad b} = \begin{bmatrix}0 \\0 \\0 \\1 \\0 \\0\end{bmatrix}}}} & \text{(11a)}\end{matrix}$

[0141] Therefore, the solution for the model parameters given to and Wcan be stated as linear least squares problem for x according to q=Ax,subject to the linear inequality constraints Gx≧b wherein q is a(m+n+p)×1 vector, A is a (m+n+p)×4 matrix, b is a 6×1 vector, and x is a4×1 vector of the model parameters. In one embodiment, A has a rank of4, i.e., full column rank and has the following singular valuedecomposition: ${A = {{U\begin{bmatrix}S \\O\end{bmatrix}}V^{T}}},$

[0142] where S is a 4×4 diagonal matrix with its diagonal entries equalto the singular values of A, U and V are orthogonal matrices ofdimensions (m+n+p)×(m+n+p) and 4×4 respectively.

[0143] The least squares solution, {circumflex over (x)}, of the leastsquares problem Q=Ax will minimize ∥Ax−Q∥². Using the singular valuedecomposition of A, ∥Ax−Q∥² can be simplified as${{{{{Ax} - Q}}^{2} = {{{{{U\begin{bmatrix}S \\O\end{bmatrix}}V^{T}x} - Q}}^{2} = {{U^{T}{U\begin{bmatrix}S \\O\end{bmatrix}}V^{T}x} - {U^{T}Q}}}}}^{2},$

[0144] since U is orthogonal, ${= {{{{\begin{bmatrix}S \\O\end{bmatrix}V^{T}x} - {U^{T}Q}}}^{2} = {{\begin{bmatrix}{{SV}^{T}x} \\0\end{bmatrix} - \begin{bmatrix}{U_{1}^{T}Q} \\{U_{2}^{T}Q}\end{bmatrix}}}^{2}}},{{and} = {{\begin{matrix}{{{SV}^{T}x} - {U_{1}^{T}Q}} \\{{- U_{2}^{T}}Q}\end{matrix}}^{2} = {{{{{SV}^{T}x} - {U_{1}^{T}Q}}}^{2} + {{U_{2}^{T}Q}}^{2}}}}$

[0145] Since ∥U₂ ^(T)Q∥² is constant, minimizing ∥Ax−Q∥² is equivalentto minimizing ∥SV^(T)x−U₁ ^(T)Q∥². In one embodiment, variable x ischanged to z, therefore,

[0146] z=SV^(T)x−U₁ ^(T)Q

V^(T)x=S⁻¹z+S⁻¹U₁ ^(T)Q. Minimizing ∥Ax−Q∥² is the same as minimizing∥z∥². Therefore, the inequality constraint becomes:

Gx≧b

G((VV ^(T))x≧b

GV(S ⁻¹ z+S ⁻¹ U ₁ ^(T) Q)≧b

GV(S ⁻¹ z+S ⁻¹ U ₁ ^(T) Q)≧b

GVS ⁻¹ z≧b−GVS ⁻¹ U ₁ ^(T) Q

[0147] Therefore, minimizing ∥Ax−Q∥² subject to the constraint GX≧b isequivalent to minimizing ∥z∥² subject to GVS⁻¹z≧b−GVS⁻¹U₁ ^(T)Q. In oneembodiment, the coefficient matrix of the above system of linear systemswill be singular value decomposed (SVD). It can be seen that thecoefficient matrix is dependent on Q(t), this means that the SVD of thecoefficient matrix has to be repeated for each tissue residue function.The coefficient matrix is N×4 where N is the number of data points ofthe tissue residue function. The first 3 columns of each coefficientmatrix is the same being only dependent on the arterial function.Therefore, SVD of the N×3 matrix can be calculated, and for each tissueresidue function, the fourth column is appended according to:$\left\lbrack {{{- k}{\sum\limits_{j = 1}^{1}\quad {{Q(j)}\Delta \quad t}}},{k{\sum\limits_{j = 1}^{2}\quad {Q(j)\Delta \quad t}}},{k{\sum\limits_{j = 1}^{3}\quad {Q(j)\Delta \quad t}}},\cdots \quad,{k{\sum\limits_{j = 1}^{N}\quad {Q(j)\Delta \quad t}}}} \right\rbrack^{T}.$

[0148] The SVD of the N×4 coeffficient matrix is then calculated from anupdate of the SVD of the basic N×3 matrix. In one embodiment, the effectof partial volume averaging (PVA) of the arterial curve on thedetermination of F and V with the deconvolution method can bedetermined. For example, letting Q(t) be the brain curve of specificmass of contrast in tissue, C_(a)(t) be the arterial curve of contrastconcentration, R(t) be the impulse residue function, F be the bloodflow, V be the blood volume, and MTT be the mean transit time. Then, ifblood flow is stationary and CT measurement is linear with respect tocontrast concentration, by principle of linear superposition thenQ(t)=FC_(a)(t)*R(t).

[0149] Also, if the arterial curve is underestimated due to partialvolume averaging from the finite resolution of CT scanners (˜8 lp/cm)then C_(a)′(t)=k.C_(a)(t), where C_(a)′(t) is the arterial curvemeasured with partial volume averaging (PVA) and k is the multiplicativefactor due to PVA, which is less than one. Therefore,${Q(t)} = {\frac{F}{k}{C_{a}^{\prime}(t)}*{{R(t)}.}}$

[0150] The deconvolution between Q(t) and C_(a)′(t) gives$\frac{F}{k}{{R(t)}.}$

[0151] The maximum height of$\frac{F}{k}{R(t)}\quad {is}\quad \frac{F}{k}$

[0152] and the area underneath$\frac{F}{k}{R(t)}\quad {is}\quad {\frac{V}{k}.}$

[0153] This explains why when a partial volume averaged arterial curveis used in the deconvolution with the brain curve, both F and V arescaled by $\frac{1}{k}{\left( {> 1} \right).}$

[0154] Note that ${MTT} = \frac{F}{V}$

[0155] is not scaled even in this case when the arterial curve ismeasured with partial volume averaging.

[0156] In one embodiment, a correction for partial volume averaging ofthe arterial curve uses the factor k such that C_(a)′(t)=k*C_(a)(t),therefore$k = {\frac{\int_{0}^{oo}{{C_{a}^{\prime}(t)}{t}}}{\int_{0}^{oo}{{C_{a}(t)}{t}}}.}$

[0157] Note that C_(a)′(t) is the measured (known) arterial curve i.e.the integral in the numerator can be calculated, while C_(a)(t) is thetrue arterial curve which is unknown. If we assume that there existregions in veins which are free of partial volume averaging, then unlikearterial curves, venous curves can be measured accurately.

[0158] For example, let h(t) be the transit time spectrum from artery toveins, therefore C_(v)(t)=C_(a)(t)*h(t). Since h(t) is a transit timespectrum, by definition ∫₀^(oo)h(t)  t = 1.

[0159] Therefore${{\int_{0}^{oo}{{C_{v}(t)}\quad {t}}} = {\int_{0}^{oo}{{C_{a}(t)}\quad {t}}}},{{{and}\quad k} = \frac{\int_{0}^{oo}{{C_{a}^{\prime}(t)}\quad {t}}}{\int_{0}^{oo}{{C_{v}(t)}\quad {t}}}},$

[0160] where k is expressed in terms of the partial volume averagedarterial curve (C_(a)′(t)) and the venous curve C_(v)(t). Both of whichare known (measured). In one embodiment, C_(a)′(t) and C_(v)(t) do notreturn to baseline during the measurement time, therefore integration toinfinity is not desirable since C_(v)(t)=C(t)*h(t), and therefore,${C_{v}(t)} = {{C_{a}^{\prime}(t)}*{\frac{h(t)}{k}.}}$

[0161] In other words, $\frac{h(t)}{k}$

[0162] can be obtained by deconvolution between C_(v)(t) and C_(a)′(t),for example if C_(a ex)′(t) is the extrapolated arterial curve obtainedby extrapolation of the trailing slope of C_(a)′(t) with an exponentialfunction. C_(a ex)′(t) returns to baseline very rapidly such that∫_(o)^(oo)C_(a, ex)^(′)(t)t

[0163] can be easily calculated. Since $\frac{h(t)}{k}$

[0164] and C_(a ex)′(t) are both known, their convolution can becalculated as${C_{v,{ex}}(t)} = {{C_{a,{ex}}^{\prime}(t)}*{\frac{h(t)}{k}.}}$

[0165] Similar to C_(a,ex)′(t),C_(v,ex)(t) also returns rapidly to zerobaseline and${{\int_{0}^{oo}{{C_{v,{ex}}(t)}{t}}} = {\frac{1}{k}{\int_{0}^{oo}{{C_{a,{ex}}^{\prime}(t)}{t}}}}},$

[0166] or${k = {\frac{\int_{0}^{oo}{{C_{a,{ex}}^{\prime}(t)}{t}}}{\int_{0}^{oo}{{C_{v,{ex}}(t)}{t}}} = \frac{\int_{o}^{T}{{C_{a,{ex}}^{\prime}(t)}{t}}}{\int_{0}^{T}{{C_{v,{ex}}(t)}{t}}}}},$

[0167] where [0,T] is the time interval within which C_(a,ex)′(t) andC_(v,ex)(t) first increase from zero baseline and then decrease back tozero baseline again. Therefore, the time shift between C_(a,ex)′(t) andC_(v,ex)(t) will not affect the calculation of k.

[0168] While the invention has been described in terms of variousspecific embodiments, those skilled in the art will recognize that theinvention can be practiced with modification within the spirit and scopeof the claims.

What is claimed is:
 1. A method for determining tissue type, said method comprising: quantitatively determining a tissue blood flow (TBF) by deconvoluting Q(t) and C_(a)(t), where Q(t) represents a curve of specific mass of contrast, and C_(a)(t) represents an arterial curve of contrast concentration; quantitatively determining a tissue blood volume (TBV) by deconvoluting Q(t) and C_(a)(t); quantitatively determining a tissue mean transit time (TMTT) by deconvoluting Q(t) and C_(a)(t); quantitatively determining a tissue capillary permeability surface area product (TPS) by deconvoluting Q(t) and C_(a)(t); and determining a tissue type based on the TBF, the TBV, the TMTT, and the TPS.
 2. A method in accordance with claim 1 wherein quantitatively determining a TBV comprises quantitatively determining a TBV for a tissue having a blood stream containing a contrast without leaking the contrast into an interstitial space of the tissue by solving a matrix equation of Q=Ah for a vector h, wherein vector h includes a plurality of elements comprising an impulse residue function at different times, Q comprises a vector including elements comprising values of a tissue residue function at different times, A comprises a matrix formed by values of the arterial curve of contrast concentration at different times, quantitatively determining a TBF comprises quantitatively determining a TBF for the tissue by solving the matrix equation of Q=Ah for the vector h, quantitatively determining a TMTT comprises quantitatively determining a TMTT for the tissue by solving the matrix equation of Q=Ah for the vector h.
 3. A method in accordance with claim 2 further comprising determining a least squares solution for the vector h under an equality constraint.
 4. A method in accordance with claim 3 wherein determining a least squares solution comprises determining a least squares solution for the vector h under a time causality constraint and a minimum transit time constraint.
 5. A method in accordance with claim 4 wherein determining a least squares solution further comprises determining a least squares solution for the vector h under a smoothness constraint, a nonnegativity constraint, and a monotonicity constraint, wherein the smoothness constraint forces h to be smoothly varying, and the monotonicity and nonnegativity constraints force h to start at a maximum and then monotonically decreases towards a zero baseline.
 6. A method in accordance with claim 1 wherein quantitatively determining a TBV comprises quantitatively determining a TBV for a tissue having a blood stream containing a contrast with leaking the contrast into an interstitial space of the tissue by solving a matrix equation of Q=Ax from the linearization of the tissue residue function for a vector x, wherein vector x includes a plurality of elements comprising TBF, TBV, TMTT, TPS and combinations thereof, Q comprises a vector including elements comprising values of a tissue residue function at different times, A comprises a matrix formed by values of the arterial curve of contrast concentration and tissue residue function at different times and combinations thereof, quantitatively determining a TBF comprises quantitatively determining a TBF for the tissue by solving the matrix equation of Q=Ax for the vector x, quantitatively determining a TMTT comprises quantitatively determining a TMTT for the tissue by solving the matrix equation of Q=Ax for the vector x, quantitatively determining a TPS comprises quantitatively determining a TPS for the tissue by solving the matrix equation of Q=Ax for the vector x.
 7. A method in accordance with claim 6 further comprising determining a least squares solution for the TBF, the TBV, the TMTT, and the TPS under a nonnegaiivity constraint that determines a TBV, a TBF, a TMTT and a TPS for the tissue where there is leakage of the contrast from the blood stream into the interstitial space.
 8. A method in accordance with claim 1 further comprising quantitatively determining a partial volume averaging scaling factor for the arterial curve of contrast concentration by: deconvoluting the measured arterial curve of contrast concentration with a venous curve of contrast concentration to determine a transit time spectrum through a tissue of interest; extrapolating the arterial curve; and convolving the extrapolated arterial curve with the transit time spectrum to generate an extrapolated venous curve, wherein the partial volume averaging scaling factor is the ratio of an area underneath the extrapolated arterial curve to an area underneath the extrapolated venous curve.
 9. A method in accordance with claim 1 further comprising scanning with at least one of a computed tomography (CT) system and Nuclear Magnetic Resonance system (NMR) to measure the arterial curve of contrast concentration and a tissue residue function, said quantitatively determining a TBV comprises quantitatively determining a TBV for a tissue having a blood stream containing a contrast without leaking the contrast into an interstitial space of the tissue by solving a matrix equation of Q=Ah for a vector h, wherein vector h includes a plurality of elements comprising an impulse residue function at different times, Q comprises a vector including elements comprising values of the tissue residue function at different times, A comprises a matrix formed by values of the arterial curve of contrast concentration at different times, said quantitatively determining a TBF comprises quantitatively determining a TBF for the tissue by solving the matrix equation of Q=Ah for the vector h, said quantitatively determining a TMTT comprises quantitatively determining a TMTT for the tissue by solving the matrix equation of Q=Ah for the vector h.
 10. A method in accordance with claim 1 further comprising scanning with at least one of a computed tomography (CT) system and Nuclear Magnetic Resonance system (NMR) to measure the arterial curve of contrast concentration and a tissue residue function, said quantitatively determining a TBV comprises quantitatively determining a TBV for a tissue having a blood stream containing a contrast with leaking the contrast into an interstitial space of the tissue by solving a matrix equation of Q=Ax from the linearization of the tissue residue function for a vector x, wherein vector x includes a plurality of elements comprising TBF, TBV, TMTT, TPS and combinations thereof, Q comprises a vector including elements comprising values of a tissue residue function at different times, A comprises a matrix formed by values of the arterial curve of contrast concentration and tissue residue function at different times and combinations thereof, said quantitatively determining a TBF comprises quantitatively determining a TBF for the tissue by solving the matrix equation of Q=Ax for the vector x, said quantitatively determining a TMTT comprises quantitatively determining a TMTT for the tissue by solving the matrix equation of Q=Ax for the vector x, said quantitatively determining a TPS comprises quantitatively determining a TPS for the tissue by solving the matrix equation of Q=Ax for the vector x.
 11. A method in accordance with claim 1 wherein the TBF is a cerebral blood flow (CBF), the TBV is a cerebral blood volume (CBV), the TMTT is a cerebral mean transit time (CMTT), the TPS is a cerebral TPS, said determining a tissue type based on the TBF, the TBV, the TMTT, and the TPS comprises determining one of a viable tissue and a non-viable tissue based on the CBF, the CBV, the CMTT, and the cerebral TBS.
 12. A system comprising at least one of a computed tomography system and a nuclear magnetic resonance system, said system configured to: quantitatively determine a tissue blood flow (TBF) by deconvoluting Q(t) and C_(a)(t), where Q(t) represents a curve of specific mass of contrast, and C_(a)(t) represents an arterial curve of contrast concentration; quantitatively determine a tissue blood volume (TBV) by deconvoluting Q(t) and C_(a)(t); quantitatively determine a tissue mean transit time (TMTT) by deconvoluting Q(t) and C_(a)(t); quantitatively determine a tissue capillary permeability surface area product (TPS) by deconvoluting Q(t) and C_(a)(t); and determine a tissue type based on the TBF, the TBV, the TMTT, and the TPS.
 13. A system according to claim 12 further configured to: quantitatively determine a TBV for a tissue having a blood stream containing a contrast without leaking the contrast into an interstitial space of the tissue by solving a matrix equation of Q=Ah for a vector h, wherein vector h includes a plurality of elements comprising an impulse residue function at different times, Q comprises a vector including elements comprising values of a tissue residue function at different times, A comprises a matrix formed by values of the arterial curve of contrast concentration at different times; quantitatively determine a TBF for the tissue by solving the matrix equation of Q=Ah for the vector h; and quantitatively determine a TMTT for the tissue by solving the matrix equation of Q=Ah for the vector h.
 14. A system according to claim 13 further configured to determine a least squares solution for the vector h under an equality constraint.
 15. A system according to claim 14 further configured to determine a least squares solution for the vector h under a time causality constraint and a minimum transit time constraint.
 16. A system according to claim 15 further configured to determine a least squares solution for the vector h under a smoothness constraint, a nonnegativity constraint, and a monotonicity constraint, wherein the smoothness constraint forces h to be smoothly varying, and the monotonicity and nonnegativity constraints force h to start at a maximum and then monotonically decreases towards a zero baseline.
 17. A system according to claim 12 further configured to: quantitatively determine a TBV for a tissue having a blood stream containing a contrast with leaking the contrast into an interstitial space of the tissue by solving a matrix equation of Q=Ax from the linearization of the tissue residue function for a vector x, wherein vector x includes a plurality of elements comprising TBF, TBV, TMTT, TPS and combinations thereof, Q comprises a vector including elements comprising values of a tissue residue function at different times, A comprises a matrix formed by values of the arterial curve of contrast concentration and tissue residue function at different times and combinations thereof; quantitatively determine a TBF for the tissue by solving the matrix equation of Q=Ax for the vector x; quantitatively determine a TMTT for the tissue by solving the matrix equation of Q=Ax for the vector x; and quantitatively determine a TPS for the tissue by solving the matrix equation of Q=Ax for the vector x.
 18. A system according to claim 17 further configured to determine a least squares solution for the TBF, the TBV, the TMTT, and the TPS under a nonnegativity constraint that determines a TBV, a TBF, a TMTT and a TPS for the tissue where there is leakage of the contrast from the blood stream into the interstitial space.
 19. A system according to claim 12 further configured to quantitatively determine a partial volume averaging scaling factor for the arterial curve of contrast concentration by: deconvoluting the measured arterial curve of contrast concentration with a venous curve of contrast concentration to determine a transit time spectrum through a tissue of interest; extrapolating the arterial curve; and convolving the extrapolated arterial curve with the transit time spectrum to generate an extrapolated venous curve, wherein the partial volume averaging scaling factor is the ratio of an area underneath the extrapolated arterial curve to an area underneath the extrapolated venous curve.
 20. A system according to claim 12, wherein the TBF is a cerebral blood flow (CBF), the TBV is a cerebral blood volume (CBV), the TMTT is a cerebral mean transit time (CMTT), the TPS is a cerebral TPS, said system further configured to determine one of a viable tissue and a non-viable tissue based on the CBF, the CBV, the CMTT, and the cerebral TBS.
 21. A computer readable medium encoded with a program executable by a computer for processing scanning data, said program configured to instruct the computer to: quantitatively determine a tissue blood flow (TBF) by deconvoluting Q(t) and C_(a)(t), where Q(t) represents a curve of specific mass of contrast, and C_(a)(t) represents an arterial curve of contrast concentration; quantitatively determine a tissue blood volume (TBV) by deconvoluting Q(t) and C_(a)(t); quantitatively determine a tissue mean transit time (TMTT) by deconvoluting Q(t) and C_(a)(t); quantitatively determine a tissue capillary permeability surface area product (TPS) by deconvoluting Q(t) and C_(a)(t); and determine a tissue type based on the TBF, the TBV, the TMTT, and the TPS.
 22. A computer readable medium in accordance with claim 21 wherein said program is further configured to: quantitatively determine a TBV for a tissue having a blood stream containing a contrast without leaking the contrast into an interstitial space of the tissue by solving a matrix equation of Q=Ah for a vector h, wherein vector h includes a plurality of elements comprising an impulse residue function at different times, Q comprises a vector including elements comprising values of a tissue residue function at different times, A comprises a matrix formed by values of the arterial curve of contrast concentration at different times; quantitatively determine a TBF for the tissue by solving the matrix equation of Q=Ah for the vector h; and quantitatively determine a TMTT for the tissue by solving the matrix equation of Q=Ah for the vector h.
 23. A computer readable medium in accordance with claim 22 wherein said program is further configured to determine a least squares solution for the vector h under an equality constraint.
 24. A computer readable medium in accordance with claim 23 wherein said program is further configured to determine a least squares solution for the vector h under a time causality constraint and a minimum transit time constraint.
 25. A computer readable medium in accordance with claim 24 wherein said program is further configured to determine a least squares solution for the vector h under a smoothness constraint, a nonnegativity constraint, and a monotonicity constraint, wherein the smoothness constraint forces h to be smoothly varying, and the monotonicity and nonnegativity constraints force h to start at a maximum and then monotonically decreases towards a zero baseline.
 26. A computer readable medium in accordance with claim 21 wherein said program is further configured to: quantitatively determine a TBV for a tissue having a blood stream containing a contrast with leaking the contrast into an interstitial space of the tissue by solving a matrix equation of Q=Ax from the linearization of the tissue residue function for a vector x, wherein vector x includes a plurality of elements comprising TBF, TBV, TMTT, TPS and combinations thereof, Q comprises a vector including elements comprising values of a tissue residue function at different times, A comprises a matrix formed by values of the arterial curve of contrast concentration and tissue residue function at different times and combinations thereof; quantitatively determine a TBF for the tissue by solving the matrix equation of Q=Ax for the vector x; quantitatively determine a TMTT for the tissue by solving the matrix equation of Q=Ax for the vector x; and quantitatively determine a TPS for the tissue by solving the matrix equation of Q=Ax for the vector x.
 27. A computer readable medium in accordance with claim 26 wherein said program is further configured to determine a least squares solution for the TBF, the TBV, the TMTT, and the TPS under a nonnegativity constraint that determines a TBV, a TBF, a TMTT and a TPS for the tissue where there is leakage of the contrast from the blood stream into the interstitial space.
 28. A computer readable medium in accordance with claim 21 wherein to quantitatively determine a partial volume averaging scaling factor for the arterial curve of contrast concentration said program is further configured to: deconvolute the measured arterial curve of contrast concentration with a venous curve of contrast concentration to determine a transit time spectrum through a tissue of interest; extrapolate the arterial curve; and convolve the extrapolated arterial curve with the transit time spectrum to generate an extrapolated venous curve, wherein the partial volume averaging scaling factor is the ratio of an area underneath the extrapolated arterial curve to an area underneath the extrapolated venous curve.
 29. A computer readable medium in accordance with claim 21 wherein the TBF is a cerebral blood flow (CBF), the TBV is a cerebral blood volume (CBV), the TMTT is a cerebral mean transit time (CMTT), the TPS is a cerebral TPS, said program is further configured to determine one of a viable tissue and a non-viable tissue based on the CBF, the CBV, the CMTT, and the cerebral TBS.
 30. A computer readable medium in accordance with claim 21 wherein said program is further configured to: scan with at least one of a computed tomography (CT) system and Nuclear Magnetic Resonance system (NMR) to measure the arterial curve of contrast concentration and a tissue residue function; quantitatively determine a TBV for a tissue having a blood stream containing a contrast without leaking the contrast into an interstitial space of the tissue by solving a matrix equation of Q=Ah for a vector h, wherein vector h includes a plurality of elements comprising an impulse residue function at different times, Q comprises a vector including elements comprising values of a tissue residue function at different times, A comprises a matrix formed by values of the arterial curve of contrast concentration at different times; quantitatively determine a TBF for the tissue by solving the matrix equation of Q=Ah for the vector h; quantitatively determine a TMTT for the tissue by solving the matrix equation of Q=Ah for the vector h; and quantitatively determine a TPS for the tissue by solving the matrix equation of Q=Ah for the vector h.
 31. A computer readable medium in accordance with claim 21 wherein said program is further configured to: scan with at least one of a computed tomography (CT) system and Nuclear Magnetic Resonance system (NMR) to measure the arterial curve of contrast concentration and a tissue residue function; quantitatively determnine a TBV for a tissue having a blood stream containing a contrast with leaking the contrast into an interstitial space of the tissue by solving a matrix equation of Q=Ax from the linearization of the tissue residue function for a vector x, wherein vector x includes a plurality of elements comprising TBF, TBV, TMTT, TPS and combinations thereof, Q comprises a vector including elements comprising values of a tissue residue function at different times, A comprises a matrix formed by values of the arterial curve of contrast concentration and tissue residue function at different times and combinations thereof; quantitatively determine a TBF for the tissue by solving the matrix equation of Q=Ax for the vector x; quantitatively determine a TMTT for the tissue by solving the matrix equation of Q=Ax for the vector x; and quantitatively determine a TPS for the tissue by solving the matrix equation of Q=Ax for the vector x. 